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INTRODUCTION 


BACKGROUND 

This report is an extension of work previously reported (Reference 1) to initial value 
problems involving stiff systems of ordinary differential equations (ODEs). As pointed out 
by Hairer and Wanner (Reference 2), the term “stiff’ is universally applied by numerical 
analysts to certain systems but does not as yet have a generally accepted precise 
mathematical definition. An early paper by Curtiss and Hirschfelder (Reference 3) supplies 
a practical definition, that stiff systems are systems of ODEs for which certain types of 
implicit numerical methods work far better tihian explicit methods. For stiff problems, 
numerical stability tends to become a bigger constraint on step-size selection than accuracy. 
As a result, methods with limited stability regions, including all explicit methods, tend to 
require a very small step size to maintain stability, which greatly increases the 
computational work and solution time, and may even cause serious roundoff errors to 
render solution infeasible. Hairer and Wanner (Reference 2) give examples of stiff systems 
that show they may consist of a single differential equation, often (but not necessarily) 
yield solutions consisting of a transient phase followed by a smooth steady-state solution, 
are often associated with physical processes that involve components of both very fast and 
very slow rates, and can result from method-of-lines discretization of partial differential 
equations. Shampine (Reference 4) points out that explicit methods may work well within 
the transient region of an otherwise stiff problem because step size is limited by accuracy, 
not stability. He concludes by noting that there is much about the theory and practice of 
solving stiff equations that is not understood. 

Since Curtiss and Hirschfelder first introduced backwards differentiation formula 
(BDF) methods (Reference 3), there has been a considerable amount of work done to 
develop software to use these methods for stiff problems. Gear (see for example Reference 
5) first developed sophisticated solvers based on BDF methods, and these are now 
tj^ically called Gear methods. EPISODE (Reference 6), LSODE (Reference 7), and 
VODE (Reference 8) are widely used codes based on Gear’s approach and include adaptive 
step-size selection and also adaptive order selection. Variable order is important because 
the BDF methods are not A-stable above second order. The entire negative real axis is 
included through sixth order, beyond which the methods are not even zero-stable, but the 
stability region becomes increasingly restricted in the area of the left half of the complex 
plane near the imaginary axis. Thus higher order methods may be used in regions where 
the Jacobian of the the derivative function has eigenvalues that are real or have small 
imaginary parts, which makes these methods effective with method-of-lines solution of 
parabolic partial differential equations. Reduction to second order is necessary for stiff 
oscillatory problems where tihe imaginary parts of the eigenvalues of the Jacobian 
dominate. Thus these methods can become very ineffective in solving hyperbolic partial 
differential equations, for example. And sixth order is considered to have a sufficiently 
restrictive stability region that it is not used in BDF codes such as LSODE (Reference 7). 
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Explicit methods of high order (8-12) have been routinely used for nonstiff problems to 
produce efficient solutions of high accuracy. Thus restriction to order two for stiff 
oscillatory problems and maximum order of five under optimal circumstances has not been 
accepted as representing true limitations, and efforts have been made to find more efficient 
stiff methods. Implicit Runge-Kutta methods were among the first to be developed. High- 
order methods with high stability have been derived, but the principal drawback has been 
the amount of work required to solve the systems of nonlinear equations for computation of 
the stages. In general this is proportional to (Ms)^, where M is ffie number of equations in 
the system and s is the number of stages, which depends on both the order and class of the 
method. For example, with Gauss methods the order is 2s, for Radau methods 2s-1, and 
for Lobatto 2s-2. Singly diagonally implicit Runge-Kutta (SDIRK) methods provide a 
significant reduction in work by enabling separate solution for each stage, a savings of a 
factor of s^. But analysis showed that with the stiff problems for which they were 
intended, the low stage order of the methods reduced the observed order of the 
convergence (see for example Reference 9). Thus singly implicit Runge-Kutta (SIRK) 
methods were developed by Burrage (Reference 10) and Butcher (Reference 11, see also 
Reference 9), which require repeated linear transformations to convert back and forth to a 
lower triangular system with a constant diagonal to be solved. The transformations are 
proportion to and so require less work for large systems. The code STRIDE 
(Reference 12) utilizes this approach, but its efficiency is seriously reduced by the amount 
of transformation work that is necessary. Butcher and Cash (Reference 13) have recently 
described an approach that reduces this transformation work somewhat. 

The class of general linear methods was formulated by Burrage and Butcher (Reference 
14, see also Reference 9) to include linear multistep methods (such as BDF methods), 
Runge-Kutta methods, and a number of other methods that had been invented over the 
years, within one theoretical formulation that would also include methods yet to be 
developed. A number of subfamilies of methods of this more general family have been 
studied recently. Blended multistep methods were developed by Skeel and Kong 
(Reference 15). SECDER (Reference 16) implements second derivative multistep methods 
studied by Enright (Reference 17). Cash first found extended BDF methods (^Reference 

18) with A-stability through fourth order and A(a)-stability through ninth order, and then 
(Reference 19) produced a much more efficient modification. Cash and Considine 
(Reference 20) have now developed codes based on these modified extended backwards 
differentiation formula (MEBDF) methods. Their code MEBDF uses variable order and 
variable step size with methods of orders 2-8. Hairer and Wanner (Reference 2) describe 
multistep Runge-Kutta methods including multistep collocation methods and multistep 
methods of Radau type in particular. Jackiewicz and Tracogna (Reference 21) and 
Tracogna (Reference 22) have described two-step Runge-Kutta methods. Butcher 
(Reference 23) proposed diagonally implicit multistage integration methods (DMSIMs). A 
class of diagonally implicit single eigenvalue methods (DIMSEMs), influenced by 
Butcher’s DMSIMs, was investigated by Enenkel (Reference 24) and Enenkel and 
Jackson (Reference 25). At this time no stiff solvers based on methods that are A-stable 
beyond second order have achieved widespread acceptance as standard recommended 
codes. Because of their high stage order and certain offier desirable design features, it is 
hoped that such a standard stiff solver may be developed based on implicit DIMSIMs. 

Butcher’s initial paper on DIMSIMs in 1992 (Reference 23) laid out the essential 
elements of a new family of general linear methods. His stated purpose was to overcome 
the glaring weaknesses of existing methods, that is, lack of A-stability for high-order linear 


6 



NAWCWPNS TP 8356 


multistep methods and low stage order and high implementation costs for A-stable implicit 
Runge-Kutta methods. These methods are diagonally implicit and hence have 
computational complexity properties similar to diagonally implicit Runge-Kutta (DIRK) 
methods, but utilize additional parameters generated through fte general linear design to 
overcome the stage order (and hence stiff order) limitations of DIRK methods. Explicit 
DIMSIMs are not technicdly “diagonally implicit,” since the diagonal is actually 0, but 
have an advantage over Runge-Kutta methods in that the order barriers for explicit Runge- 
Kutta methods do not apply, and p-stage methods of order p do indeed exist for aU positive 
integers p. The original concept called for stage order q to equal the number of internal 
stages s, the number of external stages r, and the order p. In a subsequent paper with 
Jackiewicz (Reference 26) the family of DIMSIMs was extended to include adjacent 
methods for which s+l=r = q, p = qorq-i-l, s = r-i-l=q, p = qorq+l, and 
s = r = q, p = q + l. A significant step toward practical utilization was taken with another 
paper by Butcher and Jackiewicz (Reference 27) that lays out techniques for error 
estimation, interpolation, and step-size changing. Jackiewicz, Vermiglio, and Zennaro 
(Reference 28) devised an alternative step-size changing strategy and showed how 
incorporation of an additional external stage could provide a satisfactory continuous 
method. In a separate paper (Reference 29) they also showed that there exist explicit 
DIMSIMs with regularity properties not possessed by explicit Runge-Kutta methods. 
Butcher, Chartier, and Jackiewicz, in an unpublished manuscript, “Nordsieck 
representation of DIMSIMs,” recently proposed an alternative representation of DIMSIMs 
with promise of simplifying analysis and implementation. Both explicit and implicit 
DIMSIMs up to the order 8 have now been found with appropriate stability properties and 
were announced by Butcher, Jackiewicz, and Mittelmann (Reference 30), extending 
techniques described in earlier papers by Butcher and Jackiewicz (References 31 and 32). 
A wide range of implementation issues was resolved and second and fifth order explicit 
DIMSIM computer codes were developed and successfully tested in my previous report 
(Reference 1). 


A QUICK OVERVIEW OF DIMSIMS 

To make this report as self-contained as is practical, a survey of some DIMSIM results 
that are important for orientation is provided. A more extensive discussion is provided in 
(Reference 1). We consider the initial value problem of dimension M: 


^ = f{t,y), y(to) = yo> ^s[ro,r]. (1) 

at 

We define matrices A, U, B and V such that A is s x s, V is r x r, U is s x r, and B is r x s. 
Let Y be composed of the s internal stages, F be composed of the s stage derivatives 

( Fj = -t- Cjh, ly j) and be composed of the r external stage values. Then if h is the 

step size, the solution advances one step through the relationships: 
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M 


y\''^ = hZbijFj + '^Vyy^;-'K 

7=1 ;=1 


( 2 ) 


The external stages are defined through Taylor expansion so that if 


= iaijh>y'-\-i) + olh"*'), (3) 

then we must have, for some constants a^, 

>!"'= ia,;AV'’(>„) + 0((l'’+'), (4) 

;=0 ' ' 


for a method of order p. The s values are chosen initially and other method parameters 
are then determined in a way to produce high stage order, that is, so that: 


Yi=y{l„_l+hCi) + o(h‘^^^), i = l,2,...,s. 


(5) 


where q is defined to be the stage order. We will restrict ourselves in this report to 
methods for which p = q = r =s. 

The conditions of Equation 2—^Equation 5 may be re-expressed in a convenient form as 
follows. Let W be the r X (p+1) matrix of tty values, and we denote the vector consisting 

of the kth column of W as a^.,. Let Z be a vector with element Zj = z^ '. Then define 
w(z) = WZ. Furthermore we may define e“ as the vector 






.^22 


e 


( 6 ) 
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Then the following theorem of Butcher (Reference 23) may be used in determining the 
coefficients of the method: 

Theorem 1: A DIMSIM (Equation 2) has order p and stage order p if and only if 


e‘^^=zAe^^ + Uwiz) + 0(zP^^), 

^ (7) 

e^w(z) = zBe^^ + Vw{z) + o(z^^^^. 

The coefficient methods are expressed in a tableau consisting of the matrices A, U, B, and 
V, arranged as follows: 


‘A U' 
B V 


( 8 ) 


DIMSMs were classed by Butcher (Reference 23) according to the form of A. A Type 3 
method has A = 0. Type 2 and Type 4 methods have a single diagonal value, with Type 4 
methods having a diagonal A and Type 2 methods a lower triangular A. Type 1 methods 
have a strictly lower triangular A. And V is chosen to be of rank 1, V = ev^, and v is 
chosen such that v^e = 1. 

It turns out that a restriction on the matrix elements of B provides assurance that order 
and stage order conditions are met, according to the following theorem derived by Butcher 
(Reference 23). 

Theorem 2. Let r = s = p, Ve = e. Then the DIMSIM 

'A /■ 

B V_ 

is of order p and stage order q = p if and only if B = B,, - AB, - VBj + VA, where 
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(9) 


and where 


<Pj{x)= n(jc-cj. 

k^j 


( 10 ) 


Note that using this theorem eliminates the elements of B as free parameters when 
deriving methods. It also has the following immediate corollary, since for any specified 
vector c and matrices A and V a construction may be completed, leaving stability as the 
principal issue in deriving new methods. 

Corollary: For each integer p > 1, DIMSIMs of order p and stage order q = p exist for 
s = r = p, U = I, where s is the number of internal stages and r is the length of the external 
stage vector. 

The stability matrix for a DIMSIM is (Reference 23) 


M{z) = V + zBiI-zA)~\ 


( 11 ) 


This is easily seen from the method (Equation 2) using the standard test problem y’ = X.y, 
y(tQ) = y^. We solve the first method equation for Y and, noting that F(Y) = XY, we 
obtain, setting z = hX., 


y = (/-zA)-‘yt'’-9. 


Then we may substitute this expression in the second equation to obtain 


= [zB{I-zAy^ + 
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Thus the region of absolute stability of a DIMSIM is the region 
A = {z: w € (j{M{z)) =» w < 1 }. 

If A includes the entire open left half plane, the method is called A-stable. We also define 
the associated stability polynomial 


p(w,z) = det(w/ - M{z )), 


( 12 ) 


A method may typically be verified to be A-stable either by using the Schur criterion (see, 
for example, Lambert (Reference 33)), or by reducing the stability polynomial to a familiar 
form associated with a Runge-Kutta method known to be A-stable. 

Nordsieck techniques are used extensively in the development and implementation of 
DIMSIMs. These use a vector of derivatives scaled with the step size and usually also with 
a factor of l/(j-l)! where j is the component number of the vector. Here we omit the exra 
factor, which can be readily provided through multiplication by a constant diagonal matrix 


S = diag(l/0!,l/l!,1/2!,...,1/p!). 


and use the term “Nordsieck vector” to refer to a closely related vector that frequently 
appears here. We define the Nordsieck vector of length p + 1 as 



(13) 


We also use the term to refer to the computed approximation to the exact Nordsieck vector. 

Two matrices are defined which are used to relate the Nordsieck vector to the internal 
and external stages of the method. Let be a vector with kth component 

Then we can find matrices B and V such that 

y(r„) = -h (14) 
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These matrices can be calculated using the following theorem, aimounced in Butcher 
(Reference 23) and proven rigorously in Butcher and Jackiewicz (Reference 27). We first 
define z as a vector of length p + 1, 



(15) 


Theorem 3 (Butcher and Jackiewicz, Reference 27): Assume that the method (Equation 
2) has order p and stage order q = porq = p- 1. Then the approximations (Equation 12) 

are correct to j if and only if 


e^z = + Vw{z) + 


(16) 


The Nordsick vector provides an output value and a convenient interpolant, as well as a 
predictor for implicit methods. To change step size we also define the diagonal matrix 


D = diag[\,5,d'^,...,8P), (17) 


h 

where 5 = 5^- —Then after using the expression for the Nordsieck vector at t„.,, that 
^-1 

we can rescale to accommodate a change of step size by using the formula 


(18) 


We now have a modified numerical process, as follows: 

(19) 
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CALCULATING INTERNAL STAGES WITH IMPLICIT DIMSIMS 

The description “diagonally implicit” comes from the form of the matrix As. If we 
restrict A to be lower triangular with a constant along the diagonal, the complexity of the 
solution of the system of nonlinear equations determining Y in Equation 2 is greatly 
reduced. If A is dense and a standard Gaussian elimination approach is used in a modified 
Newton method, the arithmetic complexity is 0((Ns)^), or O(N^s^), where s is the number 
of stages of the method and N is the dimension of the initial value problem. Simply 
requiring A to be lower triangular separates the system of ns simultaneous nonlinear 
equations into s systems of N simultaneous equations to be solved in sequence, each with 
arithmetic complexity 0(N^) for total complexity O(sN^), a reduction by a factor of O(s^). 

The advantage of having a single value along the diagonal may be seen from a closer 
examination of Ae solution process. The Newton iteration to solve the nonlinear system 
g(x) = 0 takes the form 


Xn+l =^n 




( 20 ) 


where J is the Jacobian of g. Of course the inverse of the Jacobian is not actually calculated 
and instead a technique such as Gaussian elimination is utilized, and the process followed 
is to calculate a correction: 


■^n+1 ~ ^n+1 • 


An LU or PLU factorization of J(xJ is called for here at each iteration, which is O(N^), 
where N is the size of the system, and this is the most time-consuming step. In practice a 
new Jacobian is evaluated and G is factored only when convergence seems too slow. In 
the case of DIMSIMs, the equation for Y takes the form 


yj=hl + +Cjh,Yj') + y^" 

^=1 


( 22 ) 


with U = I as is usually the case. Then we are solving an equation of the form gj(Yj) = 0, 
where gj takes the form: 


k=l 


(23) 
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All the dependency on Yj is contained in the first two terms and so only these terms affect 
the calculation of Ae Jacobian. Now if all diagonal elements a^ are the same the jacobians I 
will vary from stage to stage only with the solution, and new jacobians and their LU 
decompositions will only have to be computed when convergence seems too slow, which 
should be relatively rare for small step size h. This can result in substantial savings. 

A similar reduction follows from utilization of a nonsingular A (perhaps dense) with a 
single eigenvalue that may be transformed into a diagonally implicit method, and this has 
been applied in the implicit Runge-Kutta code STRIDE developed by Burrage, Butcher and 
Chipman (Reference 12). Because of the work required for linear transformations, this is 
to be avoided if possible due to the number of matrix multiplications that become 
necessary. Of course the low stage order of SDIRK methods, which drastically reduces 
the observed order of the method for stiff equations to second order, made this itemative 
approach attractive for use in STRIDE. 


INTERPOLATION AND ERROR ESTIMATION 

An efficient adaptive solver requires error estimation to enable step-size control and 
interpolation to free the choice of step size from the selection of output points. A brief 
summary of results developed more extensively and tested in (Reference 1) is provided 
here as an aid to the reader. 

A family of DIMSIM interpolants was first described in Butcher and Jackiewicz 
(Reference 27), and, for those DIMSIMs for which they exist, they provide continuous 
interpolation of maximal order with no derivative function evaluation or system solution 
cost. These interpolants are of the form 


v{tn-i + BK) = 


(24) 


Here we define ^o(^) = [^oi(^)’/5o2(^X-.^o^(0)] and ro(^) = [roi(^>X7o2(6>),-,7or(0)]^ 
where the components are polynomials of degree p (or lower if certain coefficients become 
set to zero). In order for compatibility with the equation for the first component of the 

Nordsieck vector, Pq{1) and 7o(l) must be equal to the first rows of B and V, 
respectively. The following theorem derived by Butcher and Jackiewicz (Reference 27) is 
used to derive these interpolants. 

Theorem 4 (Butcher and Jackiewicz): If a DIMSIM has order p and stage order q = p or 
q = p - 1, then 7] approximates y with uniform order p if and only if 

zl5o{d)e^^ + 'Yq{6)w{z) = e^+ 0 e (0,1], (25) 


and w(z) is as defined above. Moreover, the interpolant 7} is continuous on the whole 
interval of integration if and only if 
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Poi0) = 0, 

roiO)WDB = pM (26) 

7o(0)WDF = ro(l). 


Very few Type 2 methods that have otherwise been found provide these interpolants, 
and it was shown in the previous report (Reference 1) that Ae Nordsieck vector may 
always be utilized for this purpose, even with methods for which the Butcher-Jackiewicz 
interpolant does not exist. 

Butcher and Jackiewicz (Reference 27) also derived error estimation results. We note 
that the following general definition for local discretization error of the external stages 
applies to all DIMSIMs. The idea is to identify what is to be called the local discretization 
error of the external stages with the term of order h*^’ in the difference between the exact 
external stages at t„ and the calculated external stages, assuming that an exact Nordsieck 
vector is used initidly at t„., and that stage and order conditions are met. 

Definition; The local discretization error ^(tj of the ith external stage y|”^of the method 
of Equation 2 at the point t„ is given by 

k=0 k=l ' ' j=lk=0 


where 


= /i E + X k = (28) 

7=1 ' ' ;=l/=0 

For the case p = q = r = s, U = I, and assuming that the local condition for stage order 
holds, we make these substitutions in Equation 27 to obtain the simplified expression 


k=0 


P P 


(29) 


7=1 *=0 


In any case, the vector of values lei(t„) we designate as the local discretization error of the 
external stages. 
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The following result derived by Butcher and Jackiewicz (Reference 27) shows that the 
local truncation error is proportional to the (p + l)st derivative of the solution vector. 

Theorem 5 The local discretization error le(t„) of the external stages of Equation 2 at the 
point t„ is given by 




where 


p+^a 




k=l 


P+l-k 

k! 


BcP 

?'■ 


(31) 


To obtain a variable step estimator, we define h„ as t„-t„., and 5 = h„/h„.,, and we seek 
vectors p = P(5) and y = y( 6) such that 




It should be noted in the following modification (Reference 1) of a theorem by Butcher 
and Jackiewicz (Reference 27) that the error formula includes the effect of rescaling and 
that error estimation is not carried out by simply using a fixed step formula with a rescaled 
external stage vector as is done in interpolation. Also note that variable step size does not 
apply before the second step. The validity of the local stage order condition is assumed. 
Define a matrix 



(33) 


and a matrix 
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1 -1 i 


0 1 -1 


0 0 0 


(P + 1)! 

(-1)^ 


(34) 


Theorem 6: Assume that the method in Equation 2 with p = q = r = s, U = I, is 
implemented in variable step size mode using Ae Nordsieck technique and that V is rank 
one, Ve = e. Then 


if )8 = fi{5) and y - y(5) satisfy the system of equations 

1) r'^WDV = 0 

2) y‘^e = 0 

3) (P'^ + ^y'^WDB)e = 0 (36) 

^W^-^ + jir^WDBGti=0, i = 2,3,...,p 

5)P^^ + -^r^WDBGip^, = v\- 


For the frequently occurring case where the first row of V is and the other rows are 0, 
the first condition simplifies to y^e = 0, eliminating it as a separate condition. 

The following theorem was derived in Reference 1 for the intial step not covered by the 
Butcher-Jackiewicz theory and also may be used to produce a suitable starting step size. 

Theorem 7 (Initial Step Error Estimate): If the solution y(x) to the problem in 
Equation 1 is sufficiently smooth and the starting vector is calculated by a method correct 
up to ©(h*^^), then the error in y,, the approximation to y(Xj) calculated using the method in 
Equation 2 is 
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lte = y{ti)-yi = 




provided vectors P and y meet the following conditions: 


\)Y^e = 0 
2)P^e + Y^ai =0 
P 

4)i^=i 

,=i p! 


+ Y^oCk = 0,k = 2,...,p 


(38) 


SOME IMPLICIT DIMSIMS 


A TYPE 2 METHOD WITH BUTCHER- 
JACKIEWICZ-TYPE INTERPOLANT 

Butcher’s second order Type 2 method described in (Reference 23), although A-stable, 
was shown (Reference 1) not to have a Butcher-Jackiewicz-type interpolant. A method 
with a Butcher-Jackiewicz-type interpolant was found, however, and this process is now 
described. 

A MATHEMATICA™ program was created which incorporated the interpolant 
conditions of Theorem 4 in the search and also tested for A stability using the Schur 
criterion (Reference 33). Examination of results led to discovery of an A-stable Type 2 
method with c = [1/2,1] and with Butcher tableau 


I 0 10 

II 5 Q 1 

16 8 ^ ^ 

li _L i 1 

64 32 4 4 

107 _23 i 3 

64 32 4 4 


For this method. 
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*1 

1 

3 ■ 

8 

16 

1 

5 

15 ’ 

16 

32 


and 


"47 

17" 


"1 

3' 

64 

32 


4 

4 

0 

1 

, v= 

0 

0 

-2 

2 


0 

0 


For a Butcher-Jackiewicz-type interpolant of the form of Equation 24, we obtain 


Po = 


36 496^ 
2 64 

170 ^ 

32 


and 


70 


5 _^ ^ 

3 3 4 

_2 S9 59^ 
.3^3 4 


Using the conditions of Theorem 6 a Butcher-Jackiewicz-type error estimate was found of 
the form 




where ^(5) = 


223(5 
3 + 78’ 


Theorem 7 yield 


1 

■ 1 ■ 


1 ■ 

II 

96 

1 

. 48. 

, and 7 = 

18 

1 

18 


For the first step, the conditions of 
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est = h^TF{Y^nYrJy^^\ 


where A =^ 


■ 1 ■ 

-2 


and Y" — 

252 


The error constant for this method is -223/768, 


which is approximately 30% larger than the error constant for the second order BDF 
method. 


L-STABLE TYPE 2 SECOND ORDER METHODS 

The condition of L-stability (see for example Reference 2) requires that the method be 
A-stable, typically determined using the Schur criterion, and also that the eigenvalues of the 
stability matrix approach 0 as z approaches infinity. For the convenient c values of [0,1] 
this condition was evaluated. Note that if we specify the method matrix A as 


A = 


A 0 

a A 


and V = ev^ where v^= [Vj,V 2 ] = [v,l-v], the free parameters are v, X, and a. The stability 
polynomial is of the form 


P 2 {z)yv^ + Pi (z) w po ( 2 ) = 0 • 


In general polynomials P 2 , p,, and p,, are all of degree 2. We note that the quadratic 
formula yields eigenvalues of the form 


-Pl(z) ±^jphz) - 4P2(^)P0^ 

2P2(z) 


Then the L-stability condition will be met if the degree of P 2 is greater than the degree ofp, 
and Pq. This can be assured by calculating these coefficients (using MATHEMATICA™) 
and then finding conditions which make the coefficients of in p, and p^ equal to 0. 
These were found to yield the equations 


—A -h 2A^ + 2v — Qv — 2Av — 0 
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and 


3/1 - 4/1^ —av + 2Av = 0. 


Solving for v and a in terms of X we obtain 


A(3/l-2) ^ -2A^+6A-3 

2A-1 ’ 3A-2 


We use these values with MATHEMATICA™ to obtain a stability pol 5 aiomial and apply the 
Schur criterion along the imaginary axis. We obtain a number of expressions which 

evaluate to 0 or positive numbers, and three conditions on the parameter X for a method to 
beA-stable. These are 


-1 + 8A- 12A^ + 16A^ 

which is satisfied for A € (.1533,3.2609), 

—1 + 4A —8A^ +6A^ >0, 


which is satisfied for X >.2872, and 


-11 + 64A - lOOA^ + 64A^ - lAA"^ > 0, 

which is satisfied for A e (.261,3.2449). The intersection of these sets is the interval 
(.2872,3.2449). The error coefficient is given by the formula 

r 5-24A + 18A^ 


This has a zero near 1.08 so we choose X= 13/12 to obtain an L-stable method with an 
error constant of 1/96. We do not attempt to use the exact zero so that the stage order 
remains equal to the order and also so that an error estimate exists. The full method, along 
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"A U' 
B V 


13 

12 

0 

1 

0 

83 

13 

0 

1 

90 

12 

811 

65 

65 

9 

480 

96 

56 

56 

1091 

1703 

65 

9 

480 

1440 

56 

56 


The matrix W for computing the starting vector is given by 


W = 


1 

1 


13 

12 

181 

180 


0 

J_ 

'U 


For computation of the Nordsieck vector and rescaling we have 


rm 

13‘ 


'65 

9 * 

480 

32 


56 

56 

0 

1 

, v= 

0 

0 

-1 

1 


0 

0 


A Butcher-Jackiewicz-type interpolant does not exist so a Nordsieck or continuous 
Nordsieck interpolant should be used. Finally, the error estimation vectors are given by 


P = 


48(1 + 5) 


^811-67345-36455^+38155^ 

- ,0 


V 


755 " 


and 


r= 


2805^(1 + 5) 


(811 + 8115 + 2805^ + 2705^,-811 -8115 +105^)^. 


For the first step error estimate we find that 
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Interestingly, the first step error constant = -7/96. 

A search of other second order A-stable DIMSMs with small error constants revealed 
the method a = 6/5, X = 3/10, v = 4/5, c, = 0, C 2 = 1 with error constant -13/300. Further 

exploration turned up the DIMSIM a = 4/3, v = 3/4, X = 1/4, Cj = 0, Cj = 1, with error 
constant of -1/48,4 times better than that of the trapezoidal rule. The small error constants 
should make these competitive even with two internal stages to evaluate. 


SECOND ORDER METHODS WITH A-STABLE 
FASAL IMPLEMENTATIONS 

First approximately same as last (FASAL) DIMSIM implementations were first studied 
in the previous report (Reference 1). For methods in which c, = 0 and Cp = 1, because of 
the high stage order of DIMSEMs, there is no sacrifice in the order of accuracy through 
using the last stage derivative of the previous step in place of the first stage derivative of Ihe 
current step, a reduction of one full stage evaluation. Because of the reduction in work 
with FASAL implementations, a study was made to determine what second order 
DIMSIMs have FASAL implementations that are L-stable or at least A-stable. These 

methods must have c, = 0, C 2 = 1, leaving only three free parameters, a, X, and v, where A 
is given by 


A = 


X 

a 


0 

X 


and V = ev\ where v = [v, 1-v]. The form of the matrix B is given by Theorem 2, and for 
a second order method with p = q = r = s = 2 this requires that B = Bq - ABj - VB 2 + VA, 
where 



“1 r 


0 

I 1 

"0 O' 

II 

2 2 
_0 2_ 

II 

-1 2_ 

II 

1 1 
-2 2. 


Then we determine the stability matrix for the FASAL implementation using the techniques 
previously developed (Reference 1). We thus first produce 


A 

"0 O' 

A 

a X 


'1 0 O' 


'0 0 1' 

A = 

a X 

, B = 

B 

, u= 

0 0 1 

,and V = 

0 

[0 M 


The 3x3 stability matrix is then given by 
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^ A/ - 1 A A 

M = zB[I-zA) U + V. 

We now form the stability polynomial 


p{w,z) = det^M - w/j. 


This gives the expression 


p(w,z) = 


+ VZ- 2A.VZ -wz + 2Awz - 2vwz + ^ 

^4Avwz + 3w^z - 4A,w^z + vw^z - 2Avw^z + 2A>v^z ^ 


2(-H-Az) 


In order for a method to be L-stable, the coefficient of the highest power of w in the 
numerator must be a higher degree polynomial than any of the coefficients of other powers 
ofw. Thus we must have 


w°: v(l-2A) = 0, 

w^ —l + 2A-2v + 4Av = 0, 

w^; 3-4A + V —2Av = 0. 


From the 0th degree condition, we learn that either v = 0 or X = 1/2. But if v = 0, the 

second equation yields A. = 1/2, while the third yields A = 3/4. So we must have A = 1/2. 
The second equation will then be automatically satisfied, but the third yields the 
contradiction 1 = 0. Thus there are no L-stable p = q = s = r= 2(0,l) methods and we 
limit our search to A-stable methods. 

In order for a FASAL implementation to be A-stable, all the roots w of p must have 
|w| < 1 for Re(z) ^0, and the roots with Iwl = 1 must be simple. As noted by Burrage 
(Reference 34), the work of Chattier concerning parallel one-block methods (Reference 35) 
may be extended to include aU general linear methods. If A has a spectrum with positive 

real part (the diagonal is typically chosen as positive), then p(M(z)) will be analytic in the 
left half plane. This enables application of the maximum modulus principle to conclude that 
in the left half plane, p(M(z)) will take its largest value on the imaginary axis. Thus we 
may restrict ourselves to the imaginary axis and let z = iy, where y is real. If all the roots w 
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have \w\ < 1 along the imaginary axis and at least one point in the left half plane also 
satisfies this condition, then the region of absolute stability will include the entire left half 
plane and the method is A-stable. 

To evaluate this condition along the imaginary axis, the Schur criterion is applied, since 
a polynomial is a Schur polynomial if all the roots have modulus less than 1 (see Reference 
33). Because of the lengthy expressions that result, the work was done using 
MATHEMATICA™. We note that we may write the stability polynomial in the form 


p(w) = + g2W^ + giw + go • 


We then form 


)|C O ^ 0 ^ ^ 

K^) = 80 ^ + 82 ^ + 83- 

We have as condition 1 that |p(0)| > |p(0)|. We then form 

Pl(w) = J-[p(0)/7(w)-p(0)p(w)]. 
w 

This is a second degree polynomial in w. If p, is a also, a Schur polynomial and condition 
1 is met, then p is a Schur polynomial. We now may write p, in the form 

Pl{w) = h2W +k[W + hQ. 

We then form 


Pl{w) = hQW^ +hiw + h2. 


We have as condition 2 that |^(0)| > |pi(0)|. We then form 


P2{w) = —[^(O)pi(w) - pi(0)pi(w)]. 

w 
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If P 2 is a Schur polynomial and conditions 1 and 2 are met, then p is a Schur polynomial. 
This is a linear equation that may be solved for w and leaves only condition 3, that 

|P2(0)|>|P2(0)|- 

Each of the three conditions is a rational function in y. We first examine the 
denominators. For condition 1 this is 4(1+^,^^)- This is clearly always positive. For 
condition 2 this is 16(l+X,^y^)^ also always positive. Finally, for condition 3 we obtain a 

denominator of 32(1+X^y^)'*. Thus all the denominators are positive and the numerators 
must be investigated. And these all have only even powers of y, so we simply determine 
requirements to ensure that these coefficients are tdl positive. We list them by Schur 
condition number and by degree in y: 

Condition 1,2nd Degree: 


(-2Z + V - 2Av)(-2A - V + IXv) > 0 


Condition 1, Oth Degree: 


4>0 

Condition 2,4th Degree: 

(—2 A — 3v + 6Av)(2A — 8A^ + 3v - 6Av + 2v^ — 8Av^ + 8A^v^) > 0 
Condition 2,2nd Degree: 

4(-l + 4A + 4A^ - 2v + 8Av - - 3v^ + 12Av^ - 12A^v^) > 0 

Condition 2, Oth Degree: 

16>0 

Condition 3,8th Degree: 

(-1 + 2A)(1 + v)(2A + V - 2lbv){2X - v + 2Av)(-2A - 3v + 6Av)^ > 0 
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Condtion 3,6th Degree: 

16(-1 + 2A)(2A^ + 3Av - 8A^v + 2v^ - 5Av^ + 2A^v^ + 3v^ - 12Av^ + 12A^v^) > 0 

Condition 3,4th Degree: 

16(-l + 2A)(l-3v)>0 
Condition 3,2nd Degree: 

0>0 

Condition 3,0th Degree: 

0>0 

We immediately note that the quantity a does not appear and hence stability is determined 
only by v and X. There are a total of 6 nontrivial inequalities to satisfy. We begin with 
Condition 3,4th Degree. Here either A > ^ and v < j or A < ^ and v > j must be true. 
We then examine Condition 3, 8th Degree. The square term must be always positive. The 
product (2A + V - 2Av)(2A - V + 2Av) is the same as Condition 1, 2nd Degree. Then 

assuming this condition to be met, we must have (-1 + 2A)(1 + v) > 0. Then either A > ^ 
and v>-l or A<^ and v<-l. But since for A<^ we cannot have simultaneously 
V < -1 and V > -j, we find that we have the necessary conditions for A-stability of A > ^ 

and V € [-1,^ . Careful analysis of the remaining inequalities shows that this is also a 
sufficient condition. Condition 1,2nd Degree may be rewritten as 

(-1 - (-1 + 2A)(1 + v))(-l + (-1 + 2A)(-1 + v)) > 0. 

We now assume the requirements on X and v to note that both factors will be negative and 
thus this condition will be met. This also ensures that Condition 3, 8th Degree will be met. 
W rewrite Condition 2,2nd Degree as 
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4[-(-l + 2A)[(-1 + 2A)(1 + v)(-l + 3v) - 4] + 2 ] > 0 


Again using our requirements on X and v, we find that -1 + 3v is nonpositive while 1 + v 
and -1 + 2A, are nonnegative. Thus (-1 + 2A)(1 + v)(-l + 3v) < 0, subtracting 4 makes this 
negative, and multiplying by the nonpositive quantity -(-1 + 2X,) makes the inequality true. 

Condition 3,6th Degree will be satisfied if 

2A^ + 3Av - 8A^v + 2v^ - 5Av^ + 2??v^ + 3v^ - 12Av^ + 12A^v^ > 0. 


This may be rewritten in a more convenient form as 


X (-i + 3v)(|(-l + 2A)(2(l + v)(i(-l + 2A)(-l + 2v)-l)-l)-i) + l >0. 


We examine this from the inside outward. Since -1 + 2v is nonpositive, the innermost 
expression must be nonpositive. Since it is multiplied by a nonnegative expression, the 
next nested expression must also be nonpositive. This too is multiplied by a nonnegative 
expression and so the next nested expression is also nonpositive. But -1 + 3v is 
nonpositve, and so the whole inequality becomes true. 

The last inequality that must be satisfied is Condition 2,4th Degree. The left-hand side 
consists of a product of two factors. The first, -2X - 3v + 6^v, may be rewritten as 


(-l + 2A)(-l + 3v)-l<0 


since -1 + 3v is nonpositive and -1 + 2X, is nonnegative. The requirement for the other 
factor may be rewritten as 


(-1 + 2A)(1 + v)(2(-l + 2A)(-1 + v) - 3) -1 < 0 


Again using the restrictions on v and X we find that the innermost expression is nonpositive 
and then is multiplied by nonnegative quantities to produce a negative, from which 1 is 
subtracted. We thus have the following theorem. 
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Theorem 8 : For any DIMSIM with p = q = r= s = 2, the FASAL implementation of the 

method will be A-stable if and only if the diagonal element X. of the A matrix is not less than 
1/2 and the element v in the first column of the V matrix is in the set [-1,1/3], 

Proof; See above. ■ 

The following result shows that for a FASAL implementation the error constant can be 
as good as that of the trapezoidal rule (- 1 / 12 ) but not better. 

Theorem 9: For any DIMSIM with p = q = r = s = 2 and with an A-stable FASAL 
implementation, the error constant C<—^, with the optimal value reached for the 

diagonal element X of the A matrix equal to 1/2. 

Proof; 

The error constant for all second order DIMSIMs with Cj = 0, C 2 = 1 is given by 
^ r 5-12X-6v + 12Av 


We rewrite this as 


c=i(-i+2;i)H+v)-TL. 


Since -1 -t- v < 0, this quantity is always negative. It reaches its smallest magnitude, -1/12, 
for X = 1 / 2 . ■ 

For these optimal methods, it is not possible to calculate an error estimate of the form 
found in Theorem 6 . 

Theorem 10; DIMSIMs with Cj = 0, C 2 = 1, and X = 1/2 do not have error estimates of 
the form 


err = 


Proof; The error estimate must satisfy the conditions of Theorem 9. We define 




m 

.few. 


and 7 = 


ri(^) 

LX2(5)J' 


Imposing the condition that Y^e = 0, we find that 


72 = — 7 i. We then apply order conditions and find that there are two conditions on b2; 
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.. -a + bya -v-¥av ^ . , v + a-av 

1)- ^—2 -= 0, which gives ^— 

5 8 ^ 


6 a + + 6 ^ 12 ^^ + 6 v- 6 av l- l • 6av-6v-6a-5^ 

2)--= 0, which gives bl =- 5 - 

125 ^ ® 65 ^ 


These cannot both be satisfied for all values of 6. ■ 

However, it is possible to chose methods that do have a FASAL implementation and 
are A-stable and that have an error estimate to have an error constant arbitrarily close to that 
of the trapizoidal rule. These methods can be of Type 4 as well. 


THIRD ORDER DIMSIMS WITH A-STABLE FASAL IMPLEMENTATIONS 

Searching has been conducted for third order Type 2 DIMSIMs with c = [0,C2,1] for 
which a FASAL implementation is A-stable. To date none have been found and it might be 
conjectured that there are no A-stable FASAL implementations of DIMSIMs of order higher 
than 2. If this conjecture is true, it constitutes an unfortunate order barrier. 


A FIFTH ORDER TYPE 2 DIMSIM 

Butcher and Jackiewicz (Reference 32) found an L-stable Type 2 DIMSIM with 
p = q= r = s = 5. Implementation parameters are here derived to enable use of this 
method as a stiff solver. The stage points are given by c = [0,1/4,172,3/4,1]. The A matrix 
is 


“ .2780538411364521 

0 

0 

0 

0 

.2204522761825798 

.2780538411364521 

0 

0 

0 

2.2948198957363657 

-.6023667080712847 

.2780538411364521 

0 

0 

5.0546209011538535 

- 1.5298762183097632 

.0971191414988231 

.2780538411364521 

0 

19.3451677801081329 

- 1.4121335130997734 

- 1.8834019985178697 

.7825339554468704 

.2780538411364521 


T ^ T 

and we set V = ev and Y = cjv , where 


V = 


-.0793854651324349 

.5543175729105773 

-1.5695895491441551 

2.3320745924436815 

-.2374171510776688 
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U = L We use Theorem 2 to calculate B. This gives 


6.044855283302179 

5.853954219943505 

6.004175007913425 

6.002703177071046 

4.481882795290198 


- 2.020000467205476 

- 1.072092372634326 

- 2.014097375842605 

- 2.556003283230891 

2.672564354868939 


.03293453364122479 

- 1.839270544389963 

.6108454298803939 

3.151551366098853 

- 1.413660973235832 


.5935789859233151 

2.410922952843391 

-.963490004887004 

- 5.493514217893924 

- 8.05815479374699 


-. 2266648512058528 ' 

-.899263047489796 

-.4051827602739023 

.4481026180673921 

.909905877341711 


We use the approach developed in the first report (Reference 1) with a form for V with the 
first row of v^ and the rest of the elements equal to 0 and then apply Theorem 2 to calculate 

B. This gives 


‘ 6.044855283302179 

- 2.020000467205476 

.03293453364122479 

.5935789859233151 

. 05138898993059922 “ 

0 

0 

0 

0 

1 

1 

16 

3 

12 

-16 

25 

"T 

44 

3 

224 

3 

152 

416 

3 

140 

3 

96 

^8 

768 

-576 

160 

256 

-1024 

1536 

-1024 

256 

and W is given by Theorem 1 as 




"1 -.2780538411364521 

0 

0 

0 

0 

1 -.2485061173190319 

-.03826346028411303 

-.006085015868847461 - 

-.0005613381279595107 

-.00003711813820580275 

1 -1.470507028801533 

.1365647564495951 

.004900562818504469 

-.001619958388073782 

-.0003656404215677 

1 -3.149917665479366 

.4066191029756902 

.02777859631520006 

-.004406329750950594 

-.00169212095991602 

1 -6.110220065073812 

.929780069812273 

.0871064932281099 

-.01678258627228048 

-.00843432100142294 


Attempts to impose the conditions of Theorem 4 to derive a Butcher-Jackiewicz-type 
interpolant unfortunately led to a contradiction and thus an interpolant of this type does not 
exist for this method. The Nordsieck and continuous Nordsieck interpolants are available, 
however. The conditions of Theorem 6 were used to derive a Butcher-Jackiewicz-type 

error estimate. Free parameters were used to impose an additional condition P^e=0 (see 
Reference 27) and to ensure that there were no poles on the positive real axis (see 
Reference 1). The estimate is of the form 


est = 


where 
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5 ^ 

'LdjSJ 


j=o 


with 


do= 6.64091327569422, 
d,= 18.61705794940406, 
d2= 17.91032428174742, 
d3= 10.11460326636336, 


and 


d4= 4.180423658325783. 


Also, we calculate 


P = 


■ 13.32250820717907 ' 
-6.014932565502038 
-10.93733645614135 
4.651916710919581 
-1.022155896455255 


and 


7 = 


■ 134.9644699887695 
-134.284141177105 
0 
0 

-.6803288116644014 
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An initial error estimate was obtained using the conditions of Theorem 7. It is of the form 


est = hp^F{Y^^^yrly^^\ 


where 


^0 = 


■ -.935098307887532 ' 
.1838804364783428 
.209786992326761 
-.04441560222950508 
-.02083333013682731 


and 


70 


.1040230658993777 ' 
0 
0 
0 

-.1040230658993777 


The error coefficient for steps after the first is v^q)^ = -.000530049899772722 and for the 
first step it is calculated to be -0.000205316. All calculations were performed using 
MATHEMATICA™. 


ALTERNATIVE ERROR ESTIMATES 


Error estimates seem to be a bit more difficult for stiff problems than for nonstiff 
problems. This is because higher order terms may dominate the error in the range of 
moderate accuracy and higher hL, where L is the relevant Lipshitz constant for the 
problem. The development of some alternative approaches for error estimation is 
described. 
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BUTCHER-JACKIEWICZ-TYPE ERROR ESTIMATES 

These estimates, of the form 


est = v^le = 


are used because the external stage error propagates in the direction of v and because this 
has successfully approximated the observed locd error. An important theorem (theorem 5) 
states that the local error is given by 




where v^(Pp is then the error constant characteristic of the method and (Pp is given by 


(Pp- 


P+1 a 

I- 




p+l-k 

k\ 


BcP 
P! ■ 


The Butcher-Jackiewicz error estimation conditions are then derived by expanding the 
difference between the calculated value for the external stage vector and Ae theoretical by 
equating Taylor expansion terms to 0 up to 0(h'’) and equating the ©(h'^') term to the local 
external stage error. For variable step size the rescaling step must be included in the 

formulation. Conditions are then obtained which provide values for the vectors P and y. 
This is the approach that was used exclusively in the first report (Reference 1), along with 
some alternative equivalent formulas. Some examples of Butcher-Jackiewicz-type error 
estimates are derived in the previous section. 


USING HIGH STAGE ORDER 

The equality of the stage order q to the order p in Butcher’s original DIMSIM design 
(Reference 23) carries with it some interesting implications that have dready been described 
in the first report (Reference 1) for FASAL implementations. The ocal stage order condition 
requires that 




34 




NAWCWPNS TP 8356 


where here y denotes the local solution of the ODE through the indicated initial point. The 
internal stage values are actually used only after application of the derivative function and 
multiplied by the step size. We note that if the local stage order condition is met, 




and so it is possible where this is advantageous to substitute stage values computed using 
other methods, other step sizes, or for other steps if they are computed for the same stage 
point without creating inaccuracies. As noted in the discussion of FASAL implementations 
in the first report (Reference 1), if the results are carried forward there may be stabihty 
implications. However, error estimation typically does not entail carrying a result forward, 
and the high stage order may thus be used for this application without generating stability 
problems. It should be noted that higher order terms (0(h^^)) will be changed. In some 
cases they will be fortuitously decreased, but in others they will be increased. But error 
estimates and other calculated quantities tend to depend only on terms that are ©(h*^'). This 
is the fundamental idea that is used below to obtain new classes of no-cost DBSISIM error 
estimators. 


COMPANION METHOD ERROR ESTIMATES 

A popular approach long used with Runge-Kutta methods for estimating error for a 
method of order p is to find companion methods of adjacent orders which differ only in the 
b vector and in the addition of one extra stage. The higher order answer is taken to closely 
approximate the solution and the error is then the difference between the solution given by 
the higher order method and the solution given by the method of interest. This difference 
should be OCh**^'), with the higher order method producing hopefully a smaller error with 
principal part OCh*^^). Sometimes no additional function evaluations are required, but then 
the lower order method will typically not then use the minimum numl^r of function 
evaluations). This is the error estimation approach used with the popular Runge-Kutta- 
Fehlberg solver RKF45 (see for example Reference 36). 

For a DIMSIM with p=q = r = s, the high stage order may be used as described 
above, along with the DIMSIM feature that finding a method of the appropriate order is just 
a matter of calculating the B matrix. Thus it is possible to find a companion method which 
may be used with no extra function evaluations for which the internal stage values are 
approximately the same (differing in ©(h*^')) and for which the error constant is 0 or at 
least different. It would seem that ideally the W matrix should be the same so that the 
external stage vector has the same meaning and may be used for either method. This is 
simply a matter of using the same c values for intern^ stage points and the same A matrix, 
which then guarantees 5ie same W matrix. The vector v is then chosen to make the error 
constant 0, if possible, and otherwise at least significantly different from the error constant 
of the original method. This is done without regard to stability since it is only used for one 

step and zero stability is guaranteed by the form of V. V and V are then determined in the 

usual way as ev and , respectively. The B matrix and the B matrix must also be 
calculated for this new method, as determined in the customary way from the foregoing 
choices. Then if the step size is fixed, and if the original method is indicated by 1 and the 
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second method by 2 , beginning with the external stage vector of the previous step we 
calculate 




yf = + vf = y((„) - + 

and 

Then the error can be estimated using the difference in the Nordsieck vector first 
components: 

yj - y,” = a(^(1 ) - +(vj - vfjyl"-'! 

= (C, -C2)hP*'/P*'\l„_,)+0{hP*^). 

Thus we obtain the estimate 

+ (v! - v[)yl''-'l). (39) 


Of course the leading coefficient ratio is exactly 1 for € 3 = 0. 

A calculation for the Type 1, second order DIMSIM first developed by Butcher 
(Reference 23) and used extensively in the previous report (Reference 1) indicates that the 
same result is obtained for fixed step size as for the companion method approach outlined 
above. The method, using c = [0,1], is 
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"00 1 O’ 

2 0 0 1 

5 1 1 i 

4 4 2 2 

3 _i i i 
-4 4 2 2 . 


The error constant C, is 1/6. We choose another c = [0,1] method of the general form 


0 0 10 

2 0 0 1 

bii bi 2 V 1-v 
b2i b22 V 1-v 


Using Butcher’s theorem we obtain a B matrix of the form 


2 — 2L 

-f(-l + v) i(-l + v) 


The error constant Cj for this method is 


vV2=l2(5“H- 


We note that for explicit methods with c, = 0 and V chosen in the customary way, the first 
row of B is the same as the first row of B. Then we may calculate the same vectors P and 

y for fixed step sizes that were defined for the Butcher-Jackiewicz-type error estimate, but 
now based on the idea of subtracting the first solution from the second. 


Similarly, 


r = 


Q 

Q-Q 
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Note that these results are independent of v and are identical to those obtained using the 
Butcher-Jackiewicz approach. 

Noting that a value of v = 5/6 produces Cj = 0, we now analyze error when simply 
subtracting the first result from the second. This now yields 




and 



Thus the two approaches agree for this simple method with fixed step size, and similar 
results have been obtained for other methods. 

If step size is not fixed, it must be noted that the step change enters in as a difference 
between two methods in the rescaling process, and this must be done separately for the two 
methods or the difference will reduce to the fixed step-size formula. Thus we have 




where 




and similarly for yf. Then 
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yi - yi = hn(Bl - ^[n-l] ^ 


= K-iWDB^F Y^'' 


•vf 


^ll + evfy[" 

V J 


T 

Since v e-\, this may finally be written as 


yi -y\ = + K-\[^iWdb^ - v[wdb^)f{y^^-^^') 


+ 




y{n-2\ 


We thus propose the following alternative heuristic error estimate. 

Companion Error Estimate: Select two DIMSIMs with the same c and A but with 
different choices of v, with error constants Cj and Cj, respectively. Then the local 
trancation error of the solution at the nth step may be estimated as 


est = 


r 


Q 


IQ-C 2 




Q-Q 




where the indices 1 and 2 are used to refer to the separate methods. 

Note that this measures the error in the computed solution at the end of a step using the 
same output internal stages from the previous step and the same external stage vector from 
two steps previous. This approach differs from the Butcher-Jackiewicz approach in the 
limit of constant step size, since for the variable step approach in the case of a step ratio of 

1 , a separate rescaling would be performed since WD&^ =B^ ^ = WDB^ and V 2 . 
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COMPARISON TESTING FOR SECOND ORDER 


A Type 2 DIMSIM with a more moderately small error constant of -5/128 was selected 
for testing. This c = [0,1] method is described by the Butcher DIMSIM tableau 


25 

24 

0 

1 

0 

311 

25 

0 

1 

324 

24 

16477 

75 

225 

17 

10368 

128 

208 

208 

22093 

11275 

225 

17 

10368 

10368 

208 

208 


We also calculate 


W = 


1 

1 


25 

"24 

649 

648 


0 

il ’ 
'24. 


we choose the standard form for the rescaling matrices: 


y= 


225 

208 

0 

0 


17 

208 

0 

0 


and thus compute 


B = 


~ 16477 
10368 

0 

-1 


175 ' 

384 

1 

1 


For the Butcher-Jackiewicz-type error estimate we obtain 




55 

64(1-h 5) 



and 
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155 


104(1 + 5) 


[1-lf 


For the first step we have 



1363 
10368 ’ 



and 




For a suitable companion method with error constant 0, we find the method with 
Butcher tableau 


25 

24 

0 

1 

0 

311 

324 

25 

24 

0 

1 

1057 

5 

15 

2 

648 

8 

13 

13 

176 

365 

15 

2 

- 81 

324 

13 

13 


The W matrix is unchanged and we select 


"15 _JS\ 
13 13 


V = 


0 


0 


0 


0 


in the customary way, and compute 


B = 


1057 

648 

0 

-1 


12 

1 

1 
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The method turns out to be A-stable. 


Finally, for a companion method with a nonzero error constant, we select the method 


25 

24 

0 

1 

0 

311 

25 

0 

1 

324 

24 

26429 

75 

225 

191 

20736 

256 

416 

416 

37661 

16475 

225 

191 

20736 

20736 

416 

416 


It again has the same W matrix, and we select 


V = 


' 225 

416 

0 

0 


191 ' 

416 

0 

0 


and compute 


B = 


' 26429 

20736 

0 

-1 


575 ' 

768 

1 

1 


The error constant is -85/256 and the method is also A-stable. 

We wish to use the Prothero-Robinson problem (Reference 37), as in the first report 
(Reference 1), to test convergence of error estimate as fixed step size is reduced, behavior 
of error estimate with varying step size, and behavior with increasing stiffness. The 
Prothero-Robinson problem displays stiffness in a single equation for an appropriate choice 
of parameters. The general form of the Prothero-Robinson equation is; 


y' = A(y - p{t)) + p'it), y(to) = yo- 


(40) 


The exact solution is 


y(0 = (yo - p(^o))exp(A(t - to)) + p{t) ■ 


(41) 
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It is interesting to note that for yQ = /?(?()) the solution is simply p(t). The sine function 

was used for p(t) and the interval of interest is t e [0,20], The value for X was -2 unless 
otherwise indicated. 

A quantity r is determined for each step, where r is defined as 


r = 



-1 


Here err represents the local error as described in Equation 1, that is, the solution of the 
initial value problem of Equation 1 but beginning at a DIMSM solution point (t„.,,y„.,) 
with exact local solution at t„ given by y(t„;t„.,,y„.,), and with calculated solution at t„ of y„. 
Thus we test 




It should be noted that it is not appropriate to restart the solution in seeking to obtain the 
expression for err. A DIMSM changes character after the first step, and this must be 
preserved in the testing process. In particular the external stage vector is not recalculated 
using a starting procedure, otherwise the separate error estimate for an intial step would be 
required. Testing will only then truly reveal the behavior of the error estimator in its use 
within a solver. 

The results of a test may be indicated both graphically and using a number of statistics. 
The quantity est is used to denote the result of the use of the error estimate provided for the 
method. A graph showing err and est together is sometimes very revealing, but it is also 
vital to look at the largest error and the solution range. And a table showing the 
percentages of error estimates yielding r less than 1, 5, 10, 25, 50, and 100% has proven 
to be of great value. Of course it is expected that test results will vary greatly with the 
tolerance used and also with the problem solved and even with the problem parameters. 
Thus any testing of implementation parameters must be considered preliminary to the more 
extensive testing required of a full solver. 

Test 1 was first described by Butcher and Jackiewicz (Reference 26). A scheme was 
devised for testing the effects of rapid step changes on error estimation accuracy. An initial 

step size h^, is chosen, along with a parameter p that is varied from 1.25 to 2. The new step 
size is calculated using a ratio given by 


(-l)^sm 

r = p 



(42) 
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Thus a cyclic pattern of step sizes is used, alternately lengthening and shrinking in 
ratios from 1/p to p until the end of the interval of interest is reached. This should be 

considered to be a very stringent test, especially for higher values of p, and serves as an 
excellent preliminary check. A quick calculation similar to the one shown above for the 
Type 1 DMSIM shows that all toee approaches lead to the same error estimate for fixed 
step size without rescaling, so fixed step size tests without separate rescaling is not repeated 
for each method. 

The results shown in Figures 1 through 4 and Tables 1 through 4 were obtained with 
X = -2. The dotted line is the estimate and the solid line is the true local error. 


j^ 1 Q-+ NSteps = 150 xIO'^ NSteps = 300 




FIGURE 1. Butcher-Jackiewicz-Type, Fixed Step-Size Error Estimate. 
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TABLE 1. Butcher-Jackiewicz-Type, Fixed Step-Size Error Estimate. 


NSteps 

% < 0.01 

% < 0.05 

%<0.10 

% < 0.25 

% < 0.50 

%< 1 

150 

0.67 

1.33 

3.33 

8.67 

28.7 

86.7 


0.33 

0.33 

0.67 

2.33 

8.0 

97.3 


0.0 

0.17 

0.50 

1.33 

97.7 

99.8 

1200 

0.33 

0.92 

1.50 

87.1 

98.8 

99.5 


^ 1 g-+ NSteps = 150 


^ 1 g-5 NSteps = 300 




FIGURE 2. Companion Method Fixed Step-Size Error Estimate. 


TABLE 2. Companion Method Fixed Step-Size Error Estimate. 


NSteps 

% < 0.01 

% < 0.05 

%<0.10 

% < 0.25 

% < 0.50 

%< 1 


0.67 

2.67 

6.67 

79.3 

94.7 

98.7 


0.33 

0.67 

2.67 

92.0 

98.3 

99.7 


0 

0.50 

1.00 

8.00 

97.0 

99.2 

1200 

0 

0.25 

0.33 

1.75 

96.9 

99.7 
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+ is local error, o is estimate 

FIGURE 3. Butcher-Jackiewicz-Type Variable Step-Size Error Estimate. 


TABLE 3. Butcher-Jackiewicz-Type Variable Step-Size Error Estimate. 


Ratio 

% < 0.01 

% < 0.05 

%<0.10 

% < 0.25 

% < 0.50 

%<1 

1.25 

1.32 

5.96 

12.6 

30.5 

46.4 

62.9 

1.50 

0 

3.29 

6.58 

15.8 

42.1 

57.2 

1.75 

0 

1.97 

5.26 

12.5 

25.0 

55.9 

2.00 

0 

3.27 

5.23 

10.5 

22.2 

55.6 
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1 Q“'^Step Ratio = 1.75 




+ is local error, o is estimate 


FIGURE 4. Companion Method Variable Step-Size Error Estimate. 


TABLE 4. Companion Method Variable Step-Size Error Estimate. 


Ratio 

% < 0.01 

% < 0.05 

%<0.1 

% < 0.25 

% < 0.5 

%< 1 

1.25 

3.31 

18.54 

29.1 

51.7 

94.0 

98.7 


2.63 

11.2 

22.4 

46.1 

72.4 

99.3 

1.75 

0.66 

6.58 

15.1 

42.1 

65.1 

99.3 


1.31 

1.31 

9.80 

34.0 

62.7 

98.0 


It may be noted that only the designation “companion method error estimates” was 
used. This is because it was found that even with rescaling and variable step sizes, the 
same identical results were obtained regardless whether the third order method or the 
second order method was used as the companion method. In these tests the companion 
method error estimates with separate rescaling produced significantly better results. 
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The effect of stiffness on error estimation was also considered (Figures 5 through 8). 
This may be seen by comparing error estimation accuracy graphs for constant step size for 

values of X of -2, -10, -1(X) and -1000. The companion error estimate was used. 

^ 1 Q-+ NSteps = 150 X 1 0”^ NSteps = 300 




FIGURE 5. Fixed Step Size, X = -2. 
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y -I n"'^ NSteps = 150 X 10"^ NSteps = 300 



FIGURE 6. Fixed Step Size, = -10. 
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NSteps = 150 NSteps = 300 



0 10 20 0 10 20 



FIGURE 7. Fixed Step Size, X = -100. 
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yr 1 q- 5 NSteps = 150 X 10"“^ NSteps = 300 






FIGURE 8. Fixed Step Size, X = -1000. 


These results seem to indicate that despite final accuracy that improves with stiffness, error 
estimation deteriorates due probably to a dominance in the error of higher order terms. In 
fact the error changes sign as might be expected with a next higher term with a higher 

derivative proportional to a power of X and thus having opposite sign. The behavior for 

A, = -10 seems to show that as the step size is reduced, the error does indeed first approach 
a point where the error in the next higher term seems to approximately cancel the error term 
of degree p + 1, and then the sign changes and the estimate and the tme local error agree 
and asymptotically become equal. 


FIFTH ORDER ERROR ESTIMATION 

The Type 2 fifth order DIMSIM described above was utilized. A companion method 
with a zero error constant has the same c vector, A matrix, and W matrix. Values for v’^ are 

T 

chosen (not uniquely) to set the error constant v (p^=0. The value 
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V = 


-1.772568719302804 

1 

1 

1 

-.2274312806971961 


•J' ^ rp 

was used in this case. V = ev and V - e\v . The B matrix for the companion method 
was then determined to be 


' 4.795023599486394 

4.60412253612772 

4.754343324097641 

4.75287149325526 

3.232051111474413 


- 1.984105740930325 

- 1.036197646359175 

- 1.978202649567453 

- 2.52010855695574 

2.70845908114409 


.7671750153188124 

- 1.105030062712375 

1.345085911557982 

3.885791847776441 

-.6794204915582444 


.3573344857080197 

2.174678452628096 

- 1.199734505102299 

- 5.72975871810922 

- 8.29439929396228 


-, 2270751049413304 " 

-.899673301225274 

-.4055930140093798 

.4476923643319147 

.909495623606234 


Finally, the matrix B is calculated as 


' 4.795023599486394 

- 1.984105740930325 

.7671750153188124 

.3573344857080197 

. 0509787361951217 ' 

0 

0 

0 

0 

1 

1 

16 

3 

12 

-16 

25 

3 

44 

3 

224 

3 

152 

416 

3 

140 

3 

96 

-448 

768 

-576 

160 

256 

-1024 

1536 

-1024 

256 


It should be noted that although the error constant is 0 for steps after the first, the error 
constant for the first step is a fairly small nonzero number, 0.0000737847. 

We now compare the error estimation resulting from using the companion method with 
the Butcher-Jackiewicz-type error estimate that we have derived. A test with constant step 
size yielded the results shown in Table 5 with increasing resolution. 


TABLE 5. Butcher-Jackiewicz-Type Error Estimate For Fixed Step Size, Order 5. 


NPts 

% < 0.01 

% < 0.05 

% <0.1 

% < 0.25 

% < 0.5 

%< 1 


0 

0 

0 

2 

3.33 

99.3 


0 

0.333 

0.333 

1.33 

4.33 

99.3 


0 

0.333 

0.333 

1.17 

80.3 

99.5 

1200 

0.167 

1.00 

1.42 

5.50 

91.8 
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End point error for these tests are shown in Table 6. 


TABLE 6. End Point Error. 


NPts 

Error 

Error/h^ 


4.60e-09 

1.09e-04 


1.67e-10 

1.27e-04 


5.68e-12 

1.38e-04 


1.65e-13 

1.28e-04 


Note that the method order clearly appears through the fixed ratio in this range of step size. 
A graph of error estimates (dotted line) versus local errors (solid line) are shown in Figure 
9. 


^ 1 g-8 NPts = 150 



^10-12 NPts = 600 




FIGURE 9. Butcher-Jackiewicz-Type Error Estimate For Fixed Step Size, Order 5. 
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TTie effect of variable step size also was tested. The problem with 600 points was 
used. The following data were obtained (Table 7). 


TABLE 7. Butcher-Jackiewicz-Type Error Estimate For Variable Step Size, Order 5. 


Ratio 

% < 0.01 

% < 0.05 

%<0.1 

% < 0.25 

% < 0.5 

%< 1 

1.25 

0.17 

0.17 

0.17 

0.17 

94.5 

99.8 

1.50 

0.33 

2.49 

5.32 

14.8 

83.6 

98.7 

1.75 

0.33 

1.99 

3.49 

12.8 

71.3 

76.4 

2.00 

0.50 

1.33 

2.65 

11.8 

66.8 



A graphical representation of the results appears in Figure 10. 


Ratio = 1.25 


1 Ratio = 1.5 






FIGURE 10. Butcher-Jackiewicz-Type Error Estimate For Variable Step Size, Order 5. 


Although the picture for changing step size looks reasonably good, the number of points 
required to achieve the error estimation accuracy desired points to difficulty in efficient 
operation of a variable step-size solver. 
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The companion method was then used for purposes of comparison. The data were 
obtained for fixed step size shown in Table 8. 


TABLE 8. Companion Method Error Estimate For Fixed Step Size, Order 5. 


Npts 

% < 0.01 

% < 0.05 

%<0.1 

% < 0.25 

% < 0.5 

%< 1 


0 

1.33 

1.33 

3.33 

7.33 

97.3 


1.33 

6.33 

13.0 

36.3 

72.0 

92.3 


0 

0.50 

1.66 

4.5 

9.33 

35.2 

1200 

0.42 

2.67 

6.42 

75.8 

94.5 

98.7 


It is interesting to note in Figure 11 that the same pattern is seen as with stiffness earlier, 
and supports the interpretation that higher order term may more easily dominate the error 
when the error constant is small. 


NPts = 150 


NPts = 300 



FIGURE 11. Companion Method Error Estimate For Fixed Step Size, Order 5. 
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Variable step size for 600 points was also explored. Note that the companion method 
produced statistically poor results for fixed step size at 600 points. The data in Table 9 and 
Figure 12 were obtained. 


TABLE 9. Companion Method Error Estimate For Variable Step Size, Order 5. 


Ratio 

% < 0.01 

% < 0.05 

% < 0.1 

% < 0.25 

% < 0.5 

%<1 

1.25 

0 

3.16 

7.15 

23.23 

41.9 

63.1 


0.17 

1.83 

4.32 

11.1 

40.7 

83.1 

1.75 

0.17 

1.50 

3.16 

13.2 

45.3 



0.33 

1.49 

3.98 

27.2 

47.1 



1 ReIio = 1.25 1 Rstio = 1.5 






FIGURE 12. Companion Method Error Estimate For Variable Step Size, Order 5. 


Astonishingly, the more the step-size changes, the more accurate the error estimation. 
This is probably to be explained in terms of the effect on step-size changing on the 
contribution to the local tmncation error on leading versus higher order terms. The testing 
here is very limited, but the companion method error estimate clearly shows promise. 
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RICHARDSON-TYPE ERROR ESTIMATION 

The local stage order condition may also be used to supply an error estimate of the 
Richardson type. So long as the stage points are evenly spaced through the interval [0,1], 
the external stage vector for two steps of length h may be rescaled to provide the extern^ 
stage input, and alternating internal stage values over two steps may be used to provide the 
approximate internal stages to yield a solution for step size 2h. We then have 




and 


A = = X'n) - C{2h)l^' ) + o(a'’+^ ). (44) 


In the second equation the notation “n,n-r’ for the internal stage vector is intended to 
indicate that every other stage point beginning with the first is used, the (p + l)st derivative 
is expanded to enable use of Ae same time point, the references n,n-l,n-2 refer to step 
numbers using the step h, and 2h is used to indicate that a quantity is used with step size 
2h. The difference then is given by 


y* - yu =( 2 '^' - 2)ca'’+V'^'(<„.2)+ o(hy*^) 


Then we may estimate the error as 


err = 



(45) 


It should be noted that two steps must be completed to obtain an estimate. Thus if an 
integration step fails because tolerance is not met, it will be necessary to repeat two steps. 
And it is vital to keep constant step size over pairs of steps. Thus this approach provides a 
less sensitive adaptive integration and can only be competitive as a result of greater 
accuracy. Finally, because the larger step size is not repeated, it is expected that this 
approach may be used for explicit or predictor-corrector methods with limited stability 
regions since only the stability of the shorter step will be a factor. 
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Testing Of Richardson Error Estimate 

Richardson estimate tests were conducted (Table 10 and Figure 13) for fixed step sizes 

with the Prothero-Robinson problem (Reference 37), and with varying X. For X values of 
smaller magnitude, results were excellent. 


TABLE 10. Richardson Error Estimate, X = -2, Order 2. 


Npts 

% < 0.01 

% < 0.05 

%<0.1 

% < 0.25 

% < 0.5 

%< 1 

150 

1.33 

8.67 

17.33 

52.67 

80.67 

92.67 


0.67 

3.33 

9.67 

81.00 

95.00 

98.33 


7.00 

91.17 

96.00 

98.67 

99.33 

99.67 

1200 

0.25 

0.75 

9.92 

99.42 

99.75 

99.75 



^ 1 g-5 NSteps = 300 



FIGURE 13. Richardson Error Estimate, X = -2, Order 2. 


As the magnitude of 1 increases, the same effects that were observed before are again seen 
in Table 11 and Figure 14. 
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TABLE 11. Richardson Error Estimate, X = -10, Order 2. 


Npts 

% < 0.01 

% < 0.05 

% < 0.1 

% < 0.25 

% < 0.5 

%<1 


0.00 

0.00 

1.33 

3.33 

5.33 



0.33 

1.00 

1.00 

3.33 

7.67 

97.00 


1.33 

10.17 

35.33 

88.83 

95.17 

98.67 

1200 

0.25 

0.75 

2.00 

87.00 

98.42 

99.50 




^10-5 NSteps = 300 



FIGURE 14. Richardson Error Estimate, X. = -10, Order 2. 


The estimate was much worse for X = -100. 

Examination of the results pointed to a different heuristic formulation as the problem 
became increasingly stiff. It appeared that the result for the double step size was four times 
better than the result for the single step size. The following formula was then tested: 
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err = 


2P-\ 


(46) 


The results as seen in Tables 12 through 16 and Figures 15 through 19 were very poor for 
small X, and smaller errors but became much better as X and the errors became larger. 


TABLE 12. Heuristic Richardson Error Estimate, X. = -10, Order 2. 


Npts 

% < 0.01 

% < 0.05 

%<0.1 

% < 0.25 

% < 0.5 

%< 1 


0.00 

0.00 

0.67 

3.33 

8.00 

99.33 


0.00 

0.00 

0.33 

0.33 

1.00 

99.67 


0.00 

0.00 

0.00 

0.00 

0.33 

98.83 

1200 

0.00 

0.00 

0.00 

0.172 

0.33 

100.00 



FIGURE 15. Heuristic Richardson Error Estimate, X, = -10, Order 2. 
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TABLE 13. Heuristic Richardson Error Estimate, A, = -100, Order 2. 


Npts 

% < 0.01 

% < 0.05 

%<0.1 

% < 0.25 

% < 0.5 

%< 1 


0.67 

4.00 

7.33 

24.00 

59.33 

85.33 


1.00 

6.00 

13.67 

60.67 

86.67 

95.00 


2.00 

13.50 

57.33 

92.50 

96.50 

98.50 

1200 

0.33 

0.33 

0.58 

0.92 

6.50 

99.75 


^1q- 5 NSteps = 150 



0 10 20 
^ 1 g-5 NSteps = 600 



^ 1 Q-+ NSteps = 300 



FIGURE 16. Heuristic Richardson Error Estimate, X. = -100, Order 2. 


TABLE 14. Heuristic Richardson Error Estimate, X = -573, Order 2. 


Npts 

% < 0.01 

% < 0.05 

% < 0.1 

% < 0.25 

% < 0.5 

%< 1 


0.67 

2.00 

5.33 

15.33 

42.00 

79.33 


0.00 

1.67 

4.00 

11.33 

54.00 

90.33 


0.17 

1.00 

2.50 

8.67 

80.33 

95.67 

1200 

0.25 

0.92 

2.00 

16.00 

95.00 

98.33 
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FIGURE 17. Heuristic Richardson Error Estimate, X = -573, Order 2. 


TABLE 15. Heuristic Richardson Error Estimate, X = -1,(X)0, Order 2. 


Npts 

% < 0.01 

% < 0.05 

% <0.1 

% < 0.25 

% < 0.5 

%< 1 

150 

0.67 

2.00 

4.67 

14.67 

40.00 

78.67 


0.33 

1.67 

4.00 

10.00 

48.00 

89.67 


0.33 

1.00 

2.00 

6.33 

67.00 

95.33 

1200 

0.17 

0.50 

1.25 

5.00 

90.92 

98.00 
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FIGURE 18. Heuristic Richardson Error Estimate X = -1,000, Order 2. 
TABLE 16. Heuristic Richardson Error Estimate, X = -10,000, Order 2. 


Npts 

% < 0.01 

% < 0.05 

%<0.1 

% < 0.25 

% < 0.5 

%< 1 


0.67 

2.67 

6.00 

14.67 

39.33 

78.00 


1.00 

1.67 

3.67 

9.33 

41.00 

88.33 


0.17 

1.00 

1.67 

5.00 

38.50 

94.33 

1200 

0.25 

0.25 

0.92 

2.50 

35.25 

97.25 
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FIGURE 19. Heuristic Richardson Error Estimate, X = -10,000, Order 2. 


We see that this formula seems to predict the error with adequate accuracy over a wide 

range of X values from -100 to -10,000. The case of X = -573 was used to eliminate the 
possibility of coincidental dependence on powers of 10. 

These numerical experiments were repeated for fifth order. Results shown in Tables 17 
through 19 and Figures 20 through 22 were excellent for X = -2, and began to deteriorate 
with increasing magnitude of X and lower accuracy. 


TABLE 17. Richardson Error Estimate, X = -2, Order 5. 


Npts 

% < 0.01 

% < 0.05 

% <0.1 

% < 0.25 

% < 0.5 

%< 1 

75 

1.33 

2.67 

4.00 

12.00 

46.67 

94.67 

150 

0.67 

2.67 

6.67 

24.67 

80.00 

96.67 


0.67 

4.00 

8.33 

57.33 

91.00 

97.67 


0.67 

4.12 

10.83 

86.17 

95.50 

98.83 
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FIGURE 20. Richardson Error Estimate, X = -2, Order 5. 
TABLE 18. Richardson Error Estimate, X = -10, Order 5. 


Npts 

% < 0.01 

% < 0.05 

% < 0.1 

% < 0.25 

% < 0.5 

%< 1 

75 

0.00 

1.33 

2.67 

8.00 

26.67 


150 

0.00 

0.67 

0.67 

2.00 

2.67 

98.00 

300 

0.00 

0.33 

1.33 

2.00 

16.67 

99.33 

600 

0.67 

0.83 

1.17 

5.17 

94.00 

99.17 
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FIGURE 21. Richardson Error Estimate, X, = -10, Order 5. 
TABLE 19. Richardson Error Estimate, X = -100, Order 5 


Npts 

% < 0.01 

% < 0.05 

% <0.1 

% < 0.25 

% < 0.5 

%< 1 

75 

0.00 

2.67 

2.67 

6.67 

13.33 

21.33 


0.00 

0.67 

0.67 

3.33 

5.33 

11.33 


0.00 

0.00 

0.00 

0.67 

0.17 

5.00 

600 

0.00 

0.00 

0.17 

0.17 

0.33 

1.17 
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1 

0.5 

0 

-0.5 

-1 


1 q- 6 NSteps = 75 



0 10 20 




FIGURE 22. Richardson Error Estimate, X = -100, Order 5. 


The results at X = -100 could only be described as very poor. Some experimentation 
resulted in a different heuristic estimate for this method that seemed to be effective with 
increasing stiffness and less accuracy, the formula 


err = 


yl-yih 

2P-1 ■ 


(47) 


The following tests shown in Tables 20 through 22 and Figures 23 through 25 were 
conducted. Fewer points were used to avoid roundoff problems at high accuracy levels. 

The value X = -5732 was chosen to eliminate coincidental dependence on powers of 10. 
Actual global convergence as shown by endpoint error ratios was sixth order for these 
problems. 
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TABLE 20. Heuristic Stiff Richardson Error Estimate, X = -1000, Order = 5. 



FIGURE 23. Heuristic Stiff Richardson Error Estimate, X = -1000, Order 5. 


TABLE 21. Heuristic Stiff Richardson Error Estimate, X = -5732, Order = 5. 


Npts 

% < 0.01 

% < 0.05 

% <0.1 

% < 0.25 

% < 0.5 

75 

1.33 

6.67 

9.33 

25.33 

41.33 


1.33 

8.00 

15.33 

37.33 

60.67 


2.33 

15.00 

29.67 

57.67 

78.33 


2.33 

12.50 

22.33 

51.17 

79.83 


%< 1 
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FIGURE 24. Heuristic Stiff Richardson Error Estimate, X = -5732, Order 5. 
TABLE 22. Heuristic Stiff Richardson Error Estimate, A, = -100, Order = 5. 


Npts 

% < 0.01 

% < 0.05 

%<0.1 

% < 0.25 

% < 0.5 

%<1 

75 

0.00 

4.00 

8.00 

21.33 

38.67 

69.33 

150 

2.00 

5.33 

10.67 

25.33 

45.33 



0.00 

1.00 

2.67 

6.67 

12.67 



2.00 

14.33 

60.50 

92.50 

96.83 
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FIGURE 25. Heuristic Stiff Richardson Error Estimate, X = -100, Order 5. 

We see that there is a range within which neither the Richardson estimate nor the 
heuristic stiff Richardson estimate is effective, but both have a range of parameters for this 
problem for which high accuracy is attained, much better than that achieved with any other 
method. Because there is no si^ificant computational cost, the Richardson approach must 
be considered as highly competitive for error estimation, especially for higher order where 
error estimation has suffered from less accuracy, and this despite the problem described 
earlier of the requirement of fixed step size across two steps. It must be said that the 
heuristic stiff estimates call for more extensive theoretical study and also testing with other 
problems with different derivative functions, including systems of equations. But the 
success for this standard problem provides encouragement *at accurate error estimation for 
DIMSIM stiff solvers is possible. 
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PREDICTOR-CORRECTOR IMPLEMENTATIONS 


ALTERNATIVE APPROACHES 

Because of the high stage order of DIMSIMs, it is possible to provide a predicted 
internal stage value for in^Ucit methods based on a Taylor series expansion using the 
Nordsieck vector, which is readily available without additional function evaluations. This 
prediction is invaluable in providing a starting point for either a modified Newton solver for 
stage values or for a function iteration approach. For a method with Cj = 0 and c^ = 1, the 
finS internal stage of the last step may instead be used as a predictor for the first stage. In 
the latter mode, which provides the basis for predictor-corrector methods, the internal stage 
formula is used to provide corrections in a fixed point iteration. 


y}$ = K 1 Oij/k-l + CjK. jj"’) + (48) 

/=1 ' ' ' ' 


Convergence testing here is based on the norm of the difference between successive stage 
value iterates. 

There are a number of possible variations on the predictor-corrector idea that may be 
pursued. There are choices with the predictor. The predictor might include a term of 
degree p + 1 by using the error estimate from the previous step along with the Nordsieck 
vector. A FASAL implementation may be used to provide the first stage if c, = 0 and 
0^=1. Since it is the derivative function of the internal stage that is used, aU types of 
iterations will terminate with an evaluation, but some fixed number m of corrector iterations 
might be used, including 0 for a PE (Predict-Evaluate) method, 1 for a PECE (Predict- 
Evaluate-Correct-Evaluate) method, and in general, PEfCE)™ for any nonnegative integer 
m. And this fixed number might be different for each stage. Methods of this type are of 
course actually explicit methods. The family of PfEC)" methods, well-known with linear 
multistep methods, is not so relevant here because it is the derivative resulting from the 
evaluation step that is the goal of the stage calculations. On the other hand, iteration to 
convergence might be pursued, yielding trae function iteration stage calculations for the 
implicit DIMSIM, but this would not be suitable for stiff problems due to difficulties in 
obtaining fixed point iteration convergence for truly stiff problems. With accurate 
prediction, the variable number of fixed point iterations required in any given calculation 
might be a very small number, however, perhaps in the range of 0 to 3. LSODE 
(Reference 7), which includes a typical implementation for predictor-corrector methods 
with iteration to convergence using linear multistep methods, allows only up to perhaps 
three iterations; and if convergence has not been achieved, the step size is shortened and the 
calculation for the step is repeated in the same way as is done when the error estimate 
indicates that required accuracy was not achieved. Note that L stability will not be an 
important criterion for the DIMSIM chosen for any of these approaches since they are 
suitable for no more than moderately stiff problems, but A-stability should be provided to 
minimize stability restrictions for the implementation. 

A very significant variant to explore that has strong potential for reducing the number of 
function evaluations required is to predict the derivative function value instead of the 
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internal stage value itself, since the Taylor series terms of the Nordsieck vector may also be 
used to calculate a derivative, and this would result in what might be described as PfCE)” 
methods. Here a derivative is predicted and the stage formula provides a corrected stage 
value. The derivative function is applied first at this point in the calculation and 
convergence is determined from the norm of the difference between successive derivative 
values scaled with the step size. It should be noted that it is only the derivative function of 
the stage values that is used and not the stage values themselves, and so this provides not 
only an application of the stage formula before any evaluations but also a convergence 
measure that directly relates to the desired result. 

If the stage values are calculated with sufficient accuracy, the DMSIM formulation for 
error estimates, Nordsieck vectors, rescaling for step-size changes, and interpolation 
should remain valid. But the stability analysis must be modified to reflect the specific form 
of the predictor and iteration process used. 


RELATIONSHIPS FOR STABILITY ANALYSIS 

Assuming a standard form for V of , we have the Nordsieck vector given by 


where explicit reference to the dependency on h is included. Some simplification is 
achieved by eliminating the troublesome reference to by using the relationship 

J -I- ev • 

We may choose the first equation and write 

where B, refers to the first row of the matrix B. Then 

y{tn-i,h) = 

Another formula is available. For if we write 
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we see that 


y[n-2] ^ y[n-l] _ j 


If the error estimate is used, the external stage vector appears in a different form: 




To use standard general linear method notation, the vector representing the external stages 
must then include the internal stages, which are also carried along for the next step. The 
stability matrix will then be (2p)x(2p). If only the derivative is predicted, the analysis is 
immediately simplified since the first component of the Nordsieck vector is not used and 

this is the only component that directly depends on 


STABILITY ANALYSIS FOR A SECOND ORDER DIMSIM 

A particular L-stable Type 2 second order DIMSM with c = [0,1] was selected for 
analysis of stability with alternative approaches to prediction and varying numbers of 
corrections. The method has Butcher tableau 


25 

24 

0 

1 

0 

311 

25 

0 

1 

324 

24 

16477 

75 

225 

17 

10368 

128 

208 

208 

22093 

11275 

225 

17 

[10368 

10368 

208 

208 


W is given by 




1 

1 


25 

24 

649 

648 


0 

13 ’ 
' 24 _ 


and the rescaling matrices are given by 
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[16477 

175 " 


'lis 

17 ■ 

10368 

384 


208 

208 

0 

1 

and V = 

0 

0 

-1 

1 


0 

0 


We note that for the rescaling matrices in this form, the second component of the Nordsieck 
vector is h times the derivative function applied to the second stage, while the third 
component is h times the difference between ^e two stages as the derivative function is 
applied. There is then no difference between using the second stage in predicting the 
derivative function of the first stage at the next step and using the Nordsieck vector. The 
error constant is -5/128. 

1) Stage 1 P. Stage 2 PCE (1 Evaluation^ 

Here the derivative of stage 2 of the previous step is used to predict the derivative at the 
first step using the FASAL property. High accuracy is expected, but the FASAL 
implementation of the DIMSIM is not A-stable. The Nordsieck vector from the previous 
step is used to obtain a prediction for the derivative of the second stage. This is then used 
with the second stage formula to provide a correction. The derivative function is then 
applied with the new stage value. We thus compute 



y[«] _ y[«-l] 


= &J2,F(rj'-'l) + a22(4''-‘l 




We note that 




SO 


y^2 ^^+>' 3 " = 

Note that single subscripts on matrices are used to denote row numbers. Then 

+ ‘‘22h(h + 
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We now use the test equation y’ = A,y and set z = hX, to yield 

^^+ 3 ^ 2 ” 

We then calculate 




We use the values calculated for the internal stages to rewrite this in terms of the previous 
external and internal stages: 










We now may assemble this into a matrix form expressing 


y^[n] 

y[«] 

h 

[n] 

y\ ^ 

[n] 




y[n 1] 
y[n-l] 

h 

[n-1] 

y\ ^ 

[n-1] 

y2 \ 


where M is the stability matrix for this predictor-corrector implementation. We identify M 
from the foregoing calculations as 


M = 


0 

2 % 2 « 22(-^1 +'^ l ) 
Z^^2^22(^l +^l) 


1 


<3212 + ^22(^2 +-^2)^ 

+^ 1 ^ + 2 ^^ 2 « 22(^2 +-% 2 ) 
^ 12 ^ 21 ^ + ^ 1 ^ + Z ^^ 22 ^ 22{^2 + ^ 2 ) 


0 

0 

V 

V 


0 

1 

Z?12Z + 1-V 
b22Z +1 - V 


For the method described above, we find that 
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0 

1 

0 

0 


_2^ 

493z 


1 


24 

162 

U 

1 

M = 

625z^ 

32954z-36975z^ 

225 

136 + 975Z 


1024 

20736 

208 

1664 


281875z^ 

3579066z-5558575z^ 

225 

11016 + 146575Z 


. 248832 

1679616 

208 

134784 


(49) 


This corresponds to the stability polynomial 




27kH;2_125i 
192 ^ 384 


(50) 


The stability region is calculated from the w roots of p. Note that 1 root is 0, leaving 3 
nonzero roots. We establish a 150x150 grid of z values in the complex plane and for each z 
calculate the root w that has the maximum modulus. A contour plot was then obtained 
corresponding to a maximum modulus of 1. The same technique will be used for each of 
the following stability region calculations. For this case the stability region was obtained as 
shown in Figure 26. 
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FIGURE 26. Stage 1 P, Stage 2 PCE. 


The real axis is included to - 0.52. 

2) Stage 1 PCE. Stage 2 PCE (2 Evaluations) 

Here the derivative of stage 2 of the previous step is used to predict the derivative at the 
first step using the FASAL property. High accuracy is expected, and one iteration with the 
stage formula is used to provide extra stability from correction and evaluation. 

The Nordsieck vector from the previous step is again used to obtain a prediction for the 
derivative of the second stage and the process is carried through as in implementation 1 
described above. We thus compute 




We again note that 
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Then 


= ha2iPiYl"^) + a22h{^ + 


We now use the test equation y’ = Xy and the value of thie first internal stage to yield, with 
z = hA,, 


i']”' = i»2l(oil4“ ''+3'!“ '^] + a22z(B2+k)y^"''^+^r'^- 


We then calculate 


yf"J = = zpyt"] + 


We use the values calculated for the internal stages to rewrite this in terms of the previous 
external and internal stages: 




v[«-l] , [«-l] 

J + y} J 





We again assemble this into a matrix form expressing 


y[«] 

h 

[n] 

yi' 

[n] 

y2. 


= M 


y[n-l] 

j^n-1] 

[/i-l] 

y\ J 

[n-1] 

.^2 ; 
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where M is the stability matrix for this predictor-corrector implementation. We identify M 
from the foregoing calculations as 


M = 


flllZ 


^2(^1 ^ 3 l)^ ^22(^2 ^32)^ ^ 21 ^ ^ 




Z^^12‘^2(4i+^ l) + 2 ^ 2022 (^ 2 +^ 2 ) ''■'■^nZ + ^12«2lZ^ b^2Z + l-V 




3 2 

*22^21^112 +^21^112 

+Z ^22^2(^2 ^ 2 ) 


V + ^>21^ + ^22^1^ ^22^ +1 - V 


(51) 


For the method described above, we find that 


0 

25z 

24 


M = 

24 

25(648z + 311z2) 

liTe 


1 

31 Iz 
324 


0 

1 


625z^ 

1024 


25(-17308z^+23325z^) 

995328 


^583200 ^ 

+856804Z 

^-303225z^y 

539136 


136 + 975Z 
1664 


281875z^ 

248832 


25(l48068z^ +3506525z^) 
80621568 


^47239200 ^ 

+93055716Z 

^-45584825z^; 

43670016 


11016-H46575Z 

134784 


This corresponds to the stability polynomial 
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p{w,z) = w'* + 


- 1 - 


1489z 

576 


65225z^ 

- w. 

124416 


54425z^ ^ 
124416 j 


+ 


^ 913z 

^576 


127325z^ ^ 
62208 ; 


w 


(52) 


The stability region was obtained as shown in Figure 27 using the same technique 
described above. 



HGURE 27. Stage 1 PCE, Stage 2 PCE. 


This includes the negative real axis to - 0.93. 

3) Sta ge 1 P. Stage 2 PECE (1 Evaluations) 

Here we use the stage value derivative from the previous step for the first stage and use 
the Nordsieck vector from the previous step to predict the value of the second stage. The 
derivative function is applied and the stage formula is used to produce a corrected value. 
Its accuracy may be tested against the predicted value in an implementation where multiple 
corrections may be used. TTie derivative function is then applied for a second time to 
produce the value to be used in subsequent calculations. We thus compute 
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corrections may be used. The derivative function is then applied for a second time to 
produce the value to be used in subsequent calculations. We thus compute 

y[”] — y[”“l] 

-h 

yW = + ha22F{y[^~^^ + 

Since 

we have 

yj""^^ + =h[^+]^ + i|3)F(yf"~^^] + 

We make a substitution to eliminate reference to . Then 

y{""^J+ 5 * 2 ””^^+ = h(^+^ + 14 )F(yt""^J j+y[”"^] - /iFiF(y[""^l). 

y|”] = /ia2iF(y]"-^l) + a22hl^(^ + ^ + - Fi)/iF(y[”-^]) + yl""^^) + y^”"^l 

We now use the test equation y’ = A,y and substitute z = hA, to yield 

y]"l = Z«2i52^”"^^ + ci22Z^(Bi + ^ + - Fi)y["-1] + a22Zyl”"^^ + ■ 

We then calculate 

y[”] = /ifiF^yt”]) + Vyt""^^ = zSyt"^ + Vyt""^]. 
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We use the values calculated for the internal stages to rewrite this in terms of the previous 
external and internal stages: 








We now may assemble this into a matrix form expressing 


ry[«]' 


y[K-l] 

■*^1 

y[n] 

h 

„[«] 

= M 

y[n-l] 

h 

[n-1] 

yi ^ 

y\ 




[n-1] 

13^2 

U 2 J 



where M is the stability matrix for this predictor-corrector implementation. We identify M 
from the foregoing calculations as 


0 


Z^bi2a22(^l 


1 ^ 22 ^ 22(^11 + 


M = 

1 0 0 

a2iZ + 

<^22(-®l2+-^2 + a22Z 1 

2 

b\2^2\^ + 

^2+ ^2^22^^+^ ^>12Z + 1-V 

2-^2 “-^2) 

^22^21^^ + + 

^^^22^22(^2 + ^2 + ^22^22^^ + ^ h.2Z +1 - V 
2^2 “-^12) 


- (53) 
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M = 



0 

1 

0 

0 

25z^ 

4976z + 13725z^ 


1 

48 

5184 

24 


(1054528Z 

-373200z^ 



625z^ 

-1029375z^) 

25(-576 + 325z^) 

136 + 975Z 

2048 

663552 

13312 

1664 


(114530112z 

-56104400z^ 

-25(-139968 


281875z^ 

-154749375z^) 

+146575z^) 

11016 + 146575Z 

497664 

53747712 

3234816 

134784 


This corresponds to the stability polynomial 


p{w,z) = w* 


+ 


- 1 + 


49z 

3M 


18775z^ ^ 
9216 ^ 


w^ + 


329z 6775 zM 2 

-+- \w 

192 4608 ) 


^ 15z 

^128 


437 ^^ 
9216 ^ 


The stability region was obtained as shown in Figure 28. 



HGURE 28. Stage 1 P, Stage 2 PECE. 


It includes the real axis to -0.71. 

4) Stage 1 P. Stage 2 PCECE r2 Evaluations! 

The same approach is used as for approach 2, but an additional correction and 
evaluation is performed for the second stage. We thus compute 

y[”] = 

A"l = + 022 ( 4 ""’' + 4 ""'’) + 4 ""''' 


We again note that 
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Then 




and 


fJ"! = +a22W^&,2,F(y|"-'))+a22A(^+j 


We now use the test equation y’ = Xy and substitute z = hX, to yield 


^2”^ = (2"21 + Z^^21«22)J2” + 2^42(4 + hY-"" + {zail +1)4” 


We then calculate 


y[”3 = + yy["~^3 = ^y[«] + 


We use the values calculated for the internal stages to rewrite this in terms of the previous 
external and internal stages: 




We now may assemble this into a matrix form expressing 


yW=^| 




(2^21+4^21022)12” ^^+ 4 a 22 (^ + ^)p^" ^^+(1 + ^^ 22 ) 4 ” 
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y[«] 

y[n] 


= M 


y[n-l] 

y[n-l] 

h 

[n-1] 

y\ ^ 

[n-l] 
^2 . 


where M is the stability matrix for this predictor-corrector implementation. We identify M 
from the foregoing calculations as 


^^^ 22(^1 +^3l) 


. M = 


^21^ + ^21^22^ 
"22(^22 +^32)^^ 


2 3 

+^ 12 ^ 21 ^ 22 ^ 


1 + 0222 


2%2^22(^21+^3l) +^3^2022(^22+^32) ^ ^22(1 +0222) + !-V 

2 3 

3 2 /n S ^ ^22^212 +^22^21^222 +^2l2 , , 

2022^22(^21+^31] 3 2/- - \ ^222(1 + 0222) +1 - V 

+2 ^22^22(^22 + -^32 j 


. (55) 


For the method described above, we find that 
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M = 


0 10 0 
625z^ 3732z + 12325z^ ^ 24 + 25z 

576 3888 24 


15625z^ 

24576 


(263632Z -(1088 

-93300z^ +7800Z 


-308125z^) 225 +8125z^) 

165888 208 13312 


(85897584Z -(264384 

-42078300z^ +3517800z 


7046875z^ -138964375z^) 225 +3664375z^) 

5971968 40310784 208 3234816 


This corresponds to the stability polynomial 


p(w,z) 


+ 


-1 + 


49z 

384 


18775z^^ 

3 . 

f 329z 6775 zM 

9216 J 

W + 

t 192 4608 J 


^ 75z 

^128 


4375z 

9216 


2 \ 

“ W . 
/ 


(56) 


It is identical to the stability polynomial for the preceding implementation (stage 1 P, stage 
2 PECE) and the stability region is of course then also identical. The stability region is then 
obtained as shown in Figure 29. 
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HGURE 29. Stage 1 P, Stage 2 PCECE. 


It includes the real axis to - 0.71. 

5) Sta ge 1 PCE. Stage 2 PtCE')^ (3 Evaluations') 

Here the derivative of stage 2 of the previous step is used to predict the derivative at the 
first step using the FASAL property. High accuracy is expected, and one iteration with the 
stage formula is used to provide extra stability from correction and evaluation. 

The Nordsieck vector from the previous step is again used to obtain a prediction for the 
derivative of the second stage. The stage formula is used to provide a corrected second 
stage value and the derivative function is applied to this result. The stage formula is used to 
provide one more correction to yield the fmal value for Yj, and the derivative function is 
applied to this result in computing the external stages. We &us compute 
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4",’ = + <I22(4”‘'' + >^2''‘‘') + 

= ha2iF^haiiF^l” + *^^ + /w22(^ + ^)i^l^”“*^J + y2” 

4"] = Aa2iF(yl"l) + a22F(l|J) + yI'-'J 
= ta2iy^ta„/|j^"-‘l) + y[”-')) 

+ ha22l{ha2iF^hai + yj""’!^ + fai22(^ + %)F(yl''-'l) + 

^ [«-l] 

+ ^2 

We now use the test equation y’ = A,y and the value of the first internal stage and set z = hX. 
to yield 

^2”^ “ (^^^21^11 + ^^^22^21<^1i)^2” 

+ (za21 + Z^^22^2l)yi"~^^ + (1 + Za22)y^2~^^ 

We then calculate 

y["^ = j+ 

We use the values calculated for the internal stages to rewrite this in terms of the previous 
external and internal stages: 

v[”~l] 1 [”~1] 

+y\ 

y["] = zB (^^^21^11 + ^^^22^21^11)^2” ^22{^ + + vy[«-l] 

+(za21 + Z^^J22«2l)yi”~^^ + (1 + ^^22)^2”'^^ 
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We again assemble this into a matrix fonn expressing 


y[«] 


= M\ 


y[n-l] 

y[n-l] 

[/2-1] 

y\ ^ 

[n-1] 

yb, ^ 


where M is the stability matrix for this predictor-corrector implementation. We identify M 
from the foregoing calculations as 


M = 


0 

^ 2(^1 + ^31 


a„z 

■'■‘^ 2(^2 ^ 2 ) 2 ^ 


1 


021^(1 + 2 ^ 22 ) 


0 

1 + Z022 


^^^12^ (^21 +'®3l) 


+Z%24(4+^2) 

+^12^22^1^112^ '*’^12^1^ (1 + 2<^2) 


^ 122(1 + 2022 ) 

V 


(57) 


^2202l0iiZ^ + ^21"112^ 

2 ^^ 22 <^ 2(^21 + ^ 3 l) + 2 ^^’ 22 ^ 2(^2 + 4 ) 
+^ 22 ^ 2 ^ 1 ^ 112 '^ 


'' + ^212 ^222(1 + 2022) 

+* 22 ^ 12^(1 + 2022 ) +1-'' 


For the method described above, we find that 
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M = 


0 


625z^ 

576 


15625z^ 

24576 


24 

25(23664z^ 

+7775z^) 

186624 


25(-527264z" 

+591600z^ 

+194375z'^) 

7962624 


25(-171795168z" 

+266811600z^ 

.4 


4313088 

(1133740800 
+2233337184z 

-1094035800z^ 


24 


31U(24 + 25z) 24 + 25Z 

7776 

(4665600 
+6854432Z 

-2425800z^ 


-(1088 

+7800Z 


-2526875z^) +8125z^) 


13312 


-(264384 

+3517800Z 


7046875z^ +87663125z^) -1139620625z^) +3664375z^) 


5971968 


1934917632 


1048080384 


3234816 


This corresponds to the stability polynomial 


p{w,z) = + 


- 1 - 




289z 40775z^ 1360625z 

576 27648 2985984 


3 ^ 


\w 


+ 


^ 287z 925z^ 3183125z^l 2 ( 625z^ 1630625z^l 

. _j- ^ -ky 


V 


576 6912 


1492992 


1024 2985984 


(58) 


The stability region was obtained as shown in Figure 30. 
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HGURE 30. Stage 1 PCE, Stage 2 PCECE. 


This includes the negative real axis to - 0.59. 

6) Sta ge 1 PCE. Stage 2 PECE (3 Evaluations') 

Here we use the stage value derivative from the previous step for the first stage to 
predict the stage value derivative for the first stage and then carry out a correction and 
evaluation using the first stage formula. The size of the change may be evaluated as a 
convergence test. We then use the Nordsieck vector from the previous step to predict the 
value of the second stage. The derivative function is applied and the stage formula is used 
to produce a corrected value. Its accuracy may be tested against the predicted value in an 
implementation where multiple corrections may be used. The derivative function is then 
applied for a second time to produce the value to be used in subsequent calculations. We 
thus compute 
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Since 

j +yy["“2] 

we have 

5*!”"^^ + 4”"^^ + 2 =h(B^+Bi + j + V V"“^^ 

We make a substitution to eliminate reference to . Then 

= h[^+R2 +j 
= )W2iF('w, iF(]^"-‘') + }{"■'>) + a22*f((^ *h*jh- «|)w(rf"-'l)+>!"■'') 

[n-1] 

+>’2 • 

We now use the test equation y’ = ?iy and substitute z = hA, to yield 

5|"' = 2^02,0, il-l"-'! + a22Z^(A + ^ + +(<»2I + <I22M"‘‘’ + 4"''' 

We then calculate 

y["] = /iFF^yt"!) + Py["“^3 = ^y[«] + 

We use the values calculated for the internal stages to rewrite this in terms of the previous 
external and internal stages: 
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>] = 




z^aiia2\Y^ + z^a22(^ + ^ 
+ z{a2i+a22)y\ ^^+^2 


+yy 


[n-l] 


We now may assemble this into a matrix fonn expressing 


y^[«] 

v[«] 

h 
y\ ^ 
yyj 


= M 


y[n-l] 

y[n-l] 

h 

[n-1] 

y{ ^ 

[n-1] 

>2 


where M is the stability matrix for this predictor-corrector implementation. We identify M 
from the foregoing calculations as 


M = 


0 


aiiz 


1 


0 


z ^ 2(^1 ^1 


2 , 

021^112 + 
< 322(^2 ^2 

2^32 “^ 2 )^^ 


{<hi +“^ 22 )^ 




buaiianz^ + + 

^^^ 12*^2 (^2 ■*■ ^2 ■*■ 

■2 ^32 “ •^ 12 ) 


^>112 + ^52(021 +^22)^^ 
+v 


*12Z 

+1-V 


(59) 


Z^^22^2(Ai ■*■ ■*■ 


3 2 

^22^1^1^21^1 

3 ^ - 

Z ^22^2 (^^2 “^22 

■2 ^32 “ •^ 2 ) 


^21^ + ^22(^21 + "22)^^ 

+v 


^22^ 

+1-V 


For the method described above, we find that 
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M = 



25z 


24 

56725z^ 

15552 


1 

1297z 

648 


-(1166400 
25(-131816z^ +1713608Z 

625z^ _+170175£)_ -1264575 z^) 

2048 1990656 1078272 

(94478400 

-25(-14316264z2 +186111432z 

281875z^ +25582975z^) -190107775z^) 

497664 161243136 87340032 


0 

1 


136 + 975Z 
1664 


11016 + 146575Z 
134784 


This corresponds to the stability polynomial 


p(w,z) = w^ + 


- 1 - 


289z 615775Z 


2 ^ 


576 248832 


w 


+ 


2B7z 53875^2 112825z^"l 


576 62208 




41472 


\w^ + 


625z^ 

1024 


36725z^ ^ 
62208 , 


w 


The stability region was obtained as shown in Figure 31. 
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FIGURE 31. Stage 1 PCE, Stage 2 PECE. 


This includes the negative real axis to - 0.48. 


STABILITY USING AN A-STABLE FASAL METHOD 


The stability regions were recalculated using the following method for which the 
FASAL implementation is A-stable to get some idea of the value of A-stability for the 
FASAL implementation in a predictor-corrector setting. The stage point vector c = [0,1], 

the error constant is —— = -0.085, and the Butcher tableau is given by 
768 


257 

0 

512 

64513 

257 

64768 

512 

387841 

1 

388608 

1536 

48575 

259 

48576 

194304 


96 




















NAWCWPNS TP 8356 


The matrix W is given by 


W = 


1 

1 


257 

512 

64511 

129536 


0 

1 

512 . 


and the rescaling matrices are given by 


387841 

388608 

385 ■ 
768 


"1 

3 

2 ‘ 

3 

0 

1 

, v = 

0 

0 

-1 

1 


0 

0 


1) Predict Stage 1. PCE Stage 2 (1 Evaluation^ 

Equation 48 applies as above, but with this method we now have 


M = 


0 

1 

0 

0 

251z 

6A161Z 

0 

1 

— 

512 

32384 



257z^ 

49643648z-64767z^ 

1 

1024-z 

786432 

49741824 

3 

1536 

66563z^ 

6292211200z + 16774653z^ 

1 

129536+ 259z 

99483648 

6292340736 

3 

194304 


This corresponds to the stability polynomial 

This gives the stability region shown in Figure 32. 


(61) 
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FIGURE 32. FASAL A-Stable, Stage 1 P, Stage 2 PCE. 


This includes the negative real axis to - 0.499. 

2) Stage 1 PCE. Stage 2 PCE (1 Evaluations') 

For the method described above, we find from Equation 50 that 
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M = 


0 

257z 

512 

257(129536z 

1 

0 

251z 

+64513z^) 

64513z 

1 

512 

33161216 

64768 



(33161216 



-257(-99157760z^ 

+99287296Z 


251z^ 

+64513z^) 

-64513z^) 

1024-z 

786432 

50935627776 

99483648 

1536 



(4194893824 



257(12617972224z^ 

+12584422400Z 


66563^2 

+16708867z^) 

+16708867z^) 

129536+ 259z 

99483648 

6443356913664 

12584681472 

194304 


This corresponds to the stability polynomial 


p{w,z) = + 


f-.- 


V 


259495Z 

129536 


2257955Z 

4521984 


2 \ 

- W 

y 


3 


+ 


129959Z 

129536 


31144673z^ ^ 
24870912 ; 


w 


4139499z" 

- w. 

16580608 


The stability region was obtained as shown in Figure 33. 
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HGURE 33. FASAL A-Stable, Stage 1 PCE, Stage 2 PCE. 

This includes the negative real axis to -1.20. 

3) Sta ge 1 P. Stage 2 PECE (2 Evaluations) 

For the method described above, we find using Equation 52 that 


K 


10 
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Af = 


0 1 0 0 

151z^ 66061312z + 66646525z" 251z 

1024 66322432 512 ^ 

(101670191104z 
-66061312z' 

151 z^ -66646525z^) 262144 - 257z^ 1024-z 

1572864 101871255552 786432 1536 

(12886448537600Z 

+17109879808z" 

66563z^ +17261449975z^) 33161216+ 66563z^ 129536+ 259z 

198967296 12886713827328 99483648 194304 


This corresponds to the stability polynomial 


p(«-,z) = w*+(-!-» 


395Q09z^ 

393216 


+ 



148289z^ \ 
196608 ) 


W 





98945z^ 

393216 



(63) 


The stability region was obtained as shown in Figure 34. 
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- 1.6 - 1.4 - 1.2 -1 - 0.8 - 0.6 - 0.4 - 0.2 0 0.2 

FIGURE 34. FASAL A-Stable, Stage 1 P, Stage 2 PECE. 

The negative real axis to -1.46 is included. 

4) Sta ge 1 P. Stage 2 PCECE (2 Evaluations') 

We evaluate for this second method, using Equation 54, 
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M = 



0 

1 

0 

0 

66049z^ 

16515328Z +16645119z^ 

0 

512 + 257Z 

262144 

16580608 

512 


(25417547776Z - 16515328z^ 



66049z^ 

-16645119z^) 

1 

524288-512z-257z" 

402653184 

25467813888 

3 

786432 


(3221612134400Z + 4277469952z^ 


(66322432+132608Z 

1710669k^ 

+431108582 

1 

+66563z^) 

50935627776 

3221678456832 

3 

99483648 


This corresponds to the stability polynomial 



The stability region was obtained as shown in Figure 35. 





FIGURE 35. FASAL A-Stable, Stage 1 P, Stage 2 PCECE 


The negative real axis to -1.46 is included. 

5) Sta ge 1 PCE. Stage 2 PtCEl^ (3 Evaluations! 

For the second method described, using Equation 56, we find that 
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HGURE 36. FASAL A-Stable, Stage 1 PCE, Stage 2 PCECE. 
The negative real axis to -1.23 is included. 

6) Stage 1 PCE. Stage 2 PECE (3 Evaluations! 

For this method, using Equation 60, we obtain 
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0 

M = 

257z 

1 

0 

257z^ 

512 

99806207z^ 

194047Z 

1 

1024 

66322432 

129536 

151z^ 

-257(-198574592z^ 
+38835Iz^) 

(66322432 

+198574592Z 

-194047z^) 

1024-z 

1572864 

101871255552 

198967296 

1536 

66563z^ 

257(25168844800z^ 
+100582909z^) 

(8389787648 

+25168844800Z 

+50258173z^) 

129536+ 259z 

198967296 

12886713827328 

25169362944 

194304 


This corresponds to the stability polynomial 


p{w,z) = w‘^ 


129453Z 

129536 


13601 117zM 3 

- \w 

9043968 J 


83z 66516139z^ 

129536 66322432 


34909158 ^^ 
397934592 j 


mP' + 


^ 251p- 
^786432 


49608967z^ ^ 
397934592 ; 


w. 


The corresponding stability region is as shown in Figure 37. 
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HGURE 37. FASAL A-Stable, Stage 1 PCE, Stage 2 PECE. 
The negative real axis is included to - 0.86. 


CONCLUSIONS 

Although only two second order DIMSIMs were studied, some directions for more 
general predictor-corrector implementations of DIMSIMs are suggested. Tlie size of the 
stability matrices (2p x 2p) and the number of alternatives make this kind of study more 
difficult for high order. 

The value of a method with an A-stable FASAL implementation was evident. The 
stability regions for the second method were significantly larger. And although a correction 
with the first stage formula improved the stability region as compared with only one 
correction iteration with the second stage, an additional iteration could better be performed 
for the second stage. But two evaluations per step seem to be required to obtain a robust 
stability region, making this approach roughly comparable to PECE implementations for 
linear multistep methods. Iteration to convergence might be tried with success to produce a 
competitive method. 

The stability regions for the first method were disappointingly small. Certainly for this 
method and perhaps in general for methods that do not have an A-stable FASAL 
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implementation, iteration to convergence is the only predictor-corrector implementation that 
might lead to a competitive solver. But the fact that the stability region for the first 
implementation was comparable to that obtained for the A-stable FASAL method is a 
surprising result that keeps this a yet somewhat open question. 

No stability advantage was seen for providing an accurate prediction of the internal 
stage value with the full Nordsieck vector as compared with predicting the internal stage 
derivative at the second stage, and in comparing the stability regions for the Stage 1-PCE, 
Stage 2-PECE and PCECE methods; the second, corresponding to predicting fce internal 
stage derivative, had the larger stability region. For iteration to convergence tfie advantage 
of concluding with an evaluation and test that are directly related to lurther computations 
would make this mode seem to be the preferred one. 

Finally, it may be possible to design special methods with large stability regions for 
predictor-corrector implementations. This issue must be deferred for future study. 


THE PROTOTYPE DIMSTIFF FAMILY OF STIFF ODE SOLVERS 


INTRODUCTION 

L-stable implicit DIMSIMs of orders 2 and 5 have been implemented to develop 
prototype versions of the computer codes DIMST1FF2 and DIMSTIFF5 to solve stiff 
systems of ordinary differential equations using adaptive step-size selection. A modified 
Newton iteration with Gaussian elimination is used to solve the nonlinear systems that arise 
in the calculation of internal stages. The user is asked to provide a Jacobian subroutine. 
Calculation of numerical approximations for Jacobians and other refinements, several of 
which are discussed below, are left for future work. 


SOME DESIGN DECISIONS FOR A STIFF SOLVER 

Some helpful recommendations made by Hairer and Wanner (Reference 2) for 
implementation of implicit Runge-Kutta methods are followed. First, because stiff 
problems are associated with derivative functions having large Lipschitz constants, errors 
due to both roundoff and especially the stage-by-stage iterative process for internal stages 
may result in very large errors. To avoid tiiese problems, we work with a quantity Z as 
follows: 





(67) 


We may then write 
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zl"> = + c,A.,z|"l + yl"-'!) 

j-1 / r 1 r n\ 

or in matrix form (for a one-dimensional system) with some shorthand notation, 

zW=/i„AF(zt"lJ . (69) 


We then note that A is diagonal and nonsingular, so 




(70) 


We use the internal stage derivative values calculated in this way for rescaling, Nordsieck 
vector computation, and error estimation instead of doing the extra work of applying the 
derivative function, which would have a large Lipschitz constant and could induce large 
errors. For a diagonally implicit method, the process of using these relationships in 
comptuting the internal stage values works out in a step-by-step manner. Derivative 
function evaluations are only performed as part of the Newton iteration for computation of 
internal stage values. For example, with a second order method and an ODE of arbitrary 
dimension, we first solve 


zf"] = -h y[” 


(71) 


We then set 

*»/(<»-! + ciKA "^+>!"■'’)= 


(72) 


and then solve 


no 
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4'"^ = K¥(fn-1 +cA,z}”4yp-l]) 

(73) 

= + ^2” '*'^ 4 '^^' 

The efficiency of an implementation depends on the frequency with which Jacobians 
are evaluated and matrix factorization is carried out. Clearly this should be done as seldom 
as possible, but also as often as necessary for convergence. This must be optimized using 
a heuristic approach based on extensive testing, which is a subject for future work. As an 
initial approach awaiting future refinement some ideas were borrowed from Hairer and 
Wanner (Reference 2) and their discussion of implementing implicit Runge-Kutta methods, 
and from Enenkel and Jackson (Reference 25) who took an dgorithm diat has been used 
successfully with the code VODE (Reference 8) and modified it to reflect differences 
between backwards differentiation formula (BDF) methods and the general linear methods 
they were using. Further modifications are made because of the differences between 
DIMSIMs and implicit Runge-Kutta methods or DIMSEMs, the latter differs most notably 
for implementation purposes in that they were designed to have all stage points coincide at 
the left grid point. Thus the end point solution for the previous step serves as a predictor 
for all internal stages, which would have similar convergence properties. But DIMSIMs 
typically have nonconfluent stage points (no two points are the same). Methods 
appropriate for waveform relaxation have these stage points distributed in some fashion 
across the interval [0,1], and for methods derived thus far these are usually equally spaced. 

Newton iteration is modified here to produce a “chord method” in which the Jacobian is 
not changed during the iteration process. Although convergence is linear rather than 
quadratic, performance of chord methods for stiff problems is considered satisfactory and 
is commonly used (see for example Reference 4). A new solution value is provided at the 
end of a fully successful step, so we re-evaluate Jacobians and carry out fresh matrix 
factorizations only at the end of successful steps, rather than to do this extra work in the 
middle of steps that may fail. At this time Jacobians are re-evaluated at a set interval, but 
this may be changed to include all times when matrix factorization is carried out. 
Factorization is currently performed when step size changes or at a maximum spacing, 
whichever occurs first. 

For methods for which c, = 0 and Cp= 1, the predictor for a 0 stage point may be either 
the final internal stage of the last step or the first component of the Nordsieck vector. A 
final choice will depend on extensive testing, but at this time, for each internal stage, the 
Nordsieck vector is used to predict the internal stage value. This is somewhat different 
from predictor-corrector methods where prediction of the derivative may result in fewer 
function evaluations. This is not expected to be effective for stiff problems because 
function iteration tends to diverge. For the DIMSTIFF codes, the Nordsieck vector is used 
to predict all stage values. 

A more accurate prediction for an internal stage value should result in fewer Newton 
iterations. If convergence has not occurred within the maximum allowable number of 
iterations, either it is time to re-evaluate the Jacobian and do a fresh matrix factorization or, 
alternatively, the step size should be shortened. At the same time, accuracy of the Taylor 
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series for the Nordsieck predictor would be expected to deteriorate away from the left grid 
point where the Nordsieck vector is provided. Also, it would seem important to quickly 
determine that a step must be repeated with shorter step size or that fresh Jacobian 
evaluation and matrix factorization is needed before too much work is done for a step, 
while after most of the stages have been successfully calculated extra effort might well be 
made to calculate the final stage or stages. Thus for DMSIMs, the number of iterations 
allowed for calculation of an internal stage should depend on how far the stage point is 
from the left grid point. At this time, for second order, the maximum allowable number of 
iterations before nonconvergence is determined is set at 1 for the left end point and 6 for the 
right end point. For fifth order the maximum allowable number of iterations is currently set 
to the same number for all steps after the first, in which only 1 is allowed. However Hairer 
and Wanner (Reference 2) indicate that greater efficiency may result through allowing more 
iterations. Finding optimal numbers for each stage is a subject for future work. 

Change of step size is an important design issue. Error estimation is used to predict the 
accuracy of a step and to evaluate its success in producing a result within tolerance. 
Careful adaptation through incremental lengthening and shortening of the step size is 
usually desirable. But the matrix to be factored is 




(74) 


where A, is the diagonal element of the DIMSIM coefficient matrix A and J is the Jacobian 
of the derivative function. This matrix is strongly dependent on the step size, and so steps 
must be varied with caution to avoid a great deal of extra work. It may not be necessary to 
re-evaluate J, but appropriate factorization may change significantly and this is a very 
expensive operation. As a result, for an initial approach, the conservative strategy is 
followed that when step size reduction is called for, the step is reduced by a factor of a 
power of 2, and step size increase is not carried out unless error estimation indicates the 
step should be at least doubled, and then a maximum factor of 2 is used. Also, step size is 
not increased for a minimum of 3 steps after Jacobian evaluation and factorization are 
carried out. 

Special consideration must be given to the strategy to be followed after a step fails. 
These failures may be due either to lack of convergence or failure of the subsequent error 
estimate to be within tolerance. Actions may include step-size reduction, fresh 
factorization, Jacobian re-evaluation, or even an error message and termination. If Newton 
convergence takes place but tolerance is not met, the approach taken here is to shrink the 
step size. If step size had just increased by a factor of 2, a return made to the original step 
size might be contemplated since factorization could be preserved and little additional work 
is needed, but this has not been implemented and requires substantial additional storage. A 
total of seven retries are allowed, and step size is shrunk by a factor of 2 each time. Fresh 
matrix factorization is required for each retry. Jacobian evaluation is carried out as well 
only if the Jacobian is not current through the end of the previous step. An error message 
is provided and processing terminates if seven retries do not suffice. If the Newton 
iteration does not converge for any stage, the step must begin again with a new step size, 
which is shrank by a factor of 4. The Jacobian is updated if it is not current. Factorization 
is performed with each retry. A maximum of 10 retries are allowed before an error 
message is provided and processing is terminated. 
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Because the penalties for step failure can be so draconian, a proactive strate^ for 
periodic updating of Jacobians and matrix factorizations is utilized. A factorization is 
carried out at least every 20 steps. If a factorization is to be performed and the Jacobian has 
not been updated for at least 40 steps, the Jacobian is recalculated first. 

Richardson error estimation is currently carried out. A logical “stiff” provides a means 
to chose either the nonstiff Richardson estimator or a heuristic Richardson-type stiff 
estimator, which differs from the nonstiff estimator only by a constant factor. Since two 
successive steps with the same step size are needed to carry out a Richardson error 
estimation, some special logic is required so that after a step-size change, no estimate is 
calculated. Also, if the second step fails, the first step becomes suspect as well since it was 
not successfully tested. Thus processing must continue at the starting point of the previous 
step. 


TESTING THE DIMSTIFF CODES 

Because of project time constraints, testing was limited to the Prothero-Robinson 

problem (Equation 40) with the differential equation parameter X assuming negative values 
of large magnitude. It was chosen especially because for this stiff problem there has been 
some success in developing heuristic Richardson-type error estimation. The starting point 
(0,1/2) was utilized. To verify that the solver had the characteristics of a stiff solver, 
comparison was made with LSODE (Reference 7). A comparison is provided here (Table 
23) first between the stiff (BDF) and nonstiff (predictor-corrector Adams-Moulton) options 

to illustrate the difference between stiff and nonstiff solvers. The parameter X. takes on a 

wide range of values, with stiffness increasing with the magnitude of X. The analytical 
Jacobian was provided. Entries are the number of steps required to integrate from 0 to 20 
with absolute tolerance of 10'’. 


TABLE 23. Stiff Versus Nonstiff Solvers on the Prothero-Robinson Problem. 


X 

AM2 

BDF2 

AM5 

BDF5 

-2 

8757 

13452 

465 

542 

-1000 

35,608 

15410 

34444 

1199 

-10,000 

386,122 

15402 

384855 

1367 


It is evident that for nonstiff problems, the nonstiff solver is more efficient. However, as 
the problem becomes increasingly stiff, the amount of work with a nonstiff solver increases 
drastically, while the work remains approximately the same for the stiff solver. 

The same problem was used to test the DIMSTIFF prototype codes. Statistics are 
provided in Table 24 on the number of steps required, number of function calls, ratio of the 
number of function calls to the number of steps times, number of internal stages per step, 
number of tolerance misses, number of convergence failures, number of times the error 
estimator was deceived (that is, error was greater than the tolerance but the step was 
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accepted), number of Jacobian evaluations, number of requests for matrix factorization, 
and relative error at the end point. Of course the Jacobian is simply the number X,, and the 
matrix is also a single number. This simply provides an indication of the operation of the 
program, which mandates that these operations be carried out with a minimum frequency as 
well as in the case of tolerance or convergence misses and step-size changes. The number 
of steps required for the LSODE BDF solver of the same order is also provided. 


TABLE 24. DIMSTIFF Solvers on the Prothero-Robinson Problem. 


Order = 2 

Order = 5 

X 

-2 

-1000 

-10,000 

-2 

-1000 

-10,000 

# steps 

10780 

24031 

21916 

356 

694 

1161 

# BDF steps 

13452 

15410 

15402 

542 

1199 

1367 

# func calls 

32353 

51250 

59671 

1970 

4215 

6605 

Func/stage 

1.50 

1.07 

1.36 

1.11 

1.21 

1.14 

# tol misses 

5 

1 

126 

0 

18 

23 

# conv fails 

0 

10 

108 

0 

2 

28 

# deceived 

1 

6 

120 

0 

9 

36 

# Jac eval 

270 

601 

548 

9 

18 


# factor 

668 

1370 

1709 

44 

97 

240 

End relerr 

6e-8 

2e-10 

2e-10 

6e-10 

2e-10 

4e-12 


We see that both DIMSTIFF2 and D1MSTIFF5 clearly behaved as stiff solvers. The 

Richardson error estimator worked very well for the nonstiff problem with X = -2, and the 
heuristic Richardson-type stiff estimators performed well enough for the other problems. 
The number of function evaluations was generally just a bit higher than the optimal 1 per 
stage, and more work may be done to improve this performance. It should be noted ^at 
additional function evaluations are required when multiple iterations are needed for 
convergence, and when steps must be repeated because iterations fail to converge or 
tolerances are missed. The low ratio indicates that prediction using the Nordsieck vector is 
working very well, but the linearity of the problem makes this an easy test for the solver. 
The fifto order solver showed the expected large improvement over the second order solver 
and actually required fewer steps than LSODE. It should be noted that the fifth order BDF 
method requires only one-fifth the work in nonlinear system solving, however. It is hoped 
that accurate error estimation and efficient step-size control will significantly reduce the 
number of DMSIM steps required. But the greatest advantage of implicit DIMSIMs will 
be seen in stiff oscillatory problems (the linearized problem has eigenvalues in the left half¬ 
plane but with large imaginary parts), where the stability regions of BDF methods with 
order above 3 become very restrictive, and in the potential for solvers of order above 5. 
The starting method worked well for second order, but more work in automatically 
selecting the h^ used in computing initial derivatives for fifth order must be done. A value 
of ho= 10'® was chosen arbitrarily and provided adequate results while the previously used 
algorithm generated an hg that was much too large. 

We may note that the problem starts out nonstiff during a transient period in which the 
true solution rapidly converges to the sine function. Here step size is determined by 
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accuracy and not stability. We would expect the nonstiff error estimator to work well 
during this region. This is followed by the bulk of the integration interval, in which the 
problem is clearly stiff and the stiff formula works well. Ideally error estimators might be 
switched back and forth, depending on stiffness. Developing stiffness detection for 
DMSIMs is a topic for future work, however. Also, it should be noted that there is a 
transition period between stiff and nonstiff regions in which both types of error estimation 
work poorly. 

These experiments clearly demonstrate the potential of using implicit DIMSIMs for 
solving stiff problems. It is evident that more work is needed on stiff error estimation 
before a generally useful stiff solver can be developed. Of course, optimization is needed 
with implementation details, such as selection of Ae hg value, determination of when to 
recompute Jacobians and refactor, and also with step-size control. We plan to pursue these 
after the larger error estimation problem is adequately resolved. Although only second and 
fifth order experimental solvers are described here, an L-stable eighth order method has 
been derived (Reference 30) and this could also be developed into a solver that would have 
exceptionally high order. The work of developing DIMSIM solvers to reach their powerful 
potential is only just beginning. 
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